arXiv: 1506.03146v3 [astro-ph.SR] 8 Jan 2017 


Draft version January 10, 2017 

Preprint typeset using ETgX style emulateapj v. 01/23/15 


MODULES FOR EXPERIMENTS IN STELLAR ASTROPHYSICS (MESA): 
BINARIES, PULSATIONS, AND EXPLOSIONS 

Bill Paxton,' Pablo Marchant,^ Josiah Schwab,^ "' Evan B. Bauer,^ Lars Bildsten,'"^ Matteo Cantiello,' 
Luc Dessart,® R. Farmer,’ H. Hu,* N. Langer,’ R.H.D. Townsend,® Dean M. Townsley,*" and F.X. Timmes’ 

Draft version January 10, 2017 


ABSTRACT 

We substantially update the capabilities of the open-source software instrument Modules for Experiments in 
Stellar Astrophysics (MESA). MESA can now simultaneously evolve an interacting pair of differentially rotating 
stars undergoing transfer and loss of mass and angular momentum, greatly enhancing the prior ability to model 
binary evolution. New MESA capabilities in fully coupled calculation of nuclear networks with hundreds of 
isotopes now allow MESA to accurately simulate advanced burning stages needed to construct supernova pro¬ 
genitor models. Implicit hydrodynamics with shocks can now be treated with MESA, enabling modeling of the 
entire massive star lifecycle, from pre-main sequence evolution to the onset of core collapse and nucleosynthe¬ 
sis from the resulting explosion. Coupling of the GYRE non-adiabatic pulsation instrument with MESA allows 
for new explorations of the instability strips for massive stars while also accelerating the astrophysical use of 
asteroseismology data. We improve treatment of mass accretion, giving more accurate and robust near-surface 
profiles. A new MESA capability to calculate weak reaction rates “on-the-fly” from input nuclear data allows 
better simulation of accretion induced collapse of massive white dwarfs and the fate of some massive stars. 
We discuss the ongoing challenge of chemical diffusion in the strongly coupled plasma regime, and exhibit 
improvements in MESA that now allow for the simulation of radiative levitation of heavy elements in hot stars. 
We close by noting that the MESA software infrastructure provides bit-for-bit consistency for all results across 
all the supported platforms, a profound enabling capability for accelerating MESA’s development. 

Keywords: stars: evolution — methods: numerical — binaries: general — stars: oscillations — nuclear reac¬ 
tions — shock waves — diffusion 
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The development of a relatively complete and quantitative 
picture of stellar evolution is one of the great drivers of as- 
trop hysics. On the observ ational side of t his impetus, the Ke¬ 
pler ( |Borucki et al.|2010 1 and the CoRoT ( Baglin et al.|2009 1 
missions continuously monitored more than 100,000 stars in 
a 100 deg^ window with a dynamic range of apparent stel¬ 
lar brightness of 10®. Highlights include the discoveries that 
nearly all y Doradus and S Scuti stars are hybrid pulsators, 
and the detect ion of solar-like oscillations in a large sample 
of red giants ( Auvergne et al. l 2009 De Ridder et al. 2009t 


Grigahcene et al. 2010[ Sedding et al. 2010[ [Christensen- 


Dalsgaard & Thompson|2011 |Chaplin & Miglio|2013| l. The 

Dark Energy Purvey is scanning 500O deg^ of tne southern 

sky in 5 optical filters every few days to discover and study 
thousands of supernovae (e.g.,|Papadopoulos et al.|2015||Yuan 
et al.|20l3] l. B uilding upon tl^ l egacy of the Palomar Tran- 
sient Factory ( |Law et al.||2009} , the intermediate Palomar 
Transient Factory conducts a fully-automated, wide-held sur¬ 


vey that systematically e xplores the transient sk y with a 90 
second to 5 day cadence ( [Vreeswijk et al.|2014| l. The forth¬ 
coming Zwicky Transient Facility will enable a survey more 
than an order of magnitude faster at the same depth as its pre¬ 
decessors. In its unique orbit the Transiting Exoplanet Survey 
Satellite will have an unobstructed view to scrutinize the light 
curves of the bright est 100,000 stars with a 1 minute cadence 
( Ricker et al.||2015 l. The Gaia mission aims to provide un- 
precedented distance and radial velocity measurements with 
the accuracies needed to reveal the evolutionary state, compo- 
sitio n, and kinematics of about one billion st ars in our Galaxy 
(e.g., Creevey et al. [2015] [Sacco et al.|2015] l. The Large Syn¬ 
optic Siirvey'THescopFwunmK^'nieennre Southern Hemi¬ 
sphere deeply in multiple optical colors every week with its 
three billion pixel digital camera, thus opening a new window 
on tr ansient objects such as interacting close binary systems, 
(e.g., Oluseyi et al.|2012) l. 

Interpreting these new observations and predicting new 
stellar phenomena propels the theoretical side, in particular 
the evolution of the community software instrument Modules 
for Experiments in Stellar Astrophysics (MESA) for research 
and education. We introduced MESA in [Paxton et al.| ( |201 1[ 
hereafter Paper I) a nd expanded its range of capabilities in 
Paxton et al.|(|2013| hereafter Paper II). This paper describes 
the major new advances for MESA modeling of binary systems, 
shock hydrodynamics, explosions of massive stars and x-ray 
bursts with large, in situ reaction networks. Moreover it de¬ 
tails the coupling of MESA with the non-adiabatic pul sation 
software instrument GYRE ( [Townsend & Teitler||2013l l. We 
also describe advances made to existing modules since Pa¬ 
per II, including improved treatments of mass accretion, weak 
reaction rates, and particle diffusion. _ 

It h as been a little more than 200 years since Herschel 
( [1802[ ) announced, after 25 years of observation, that certain 
pairs of stars displayed evidence of orbital motion around 
their common center of mass. Binary systems allow the 
masses of their component stars to be directly determined, 
which in turn allows stellar radii to be indirectly estimated. 
This allows the calibration of an empirical mass-luminosity 
relationship from which the masses of single stars can be esti¬ 


mated (Torres et al.|2010 i. Recent surveys such as Raghavaii] 


et al. (2010[l suggest 30% to 50% of solar-like systems in the 


Galactic disk are composed of binaries, where the binary frac- 

r 


tion is higher for m ore massive stars ( Sana et al.|2012 |Kob- 
ulnicky et al.|2014] l. As argued by de Mink et al.| ( ^13[ l, the 
most rapidly rotating massive stars are expected in binary sys¬ 
tems as a consequence of accretion-induced spin-up. Differ¬ 
ential rotationhasam^orin^ctontheeyohition of massive 


stars ( Maeder & Meynet|2000j[Heger et al.[[2000 2005| l and 
for single stars the corresponding physics has been included 
into MESA as described in Paper II. On the other hand very 
few works that include the ph ysics of differential rotation in 
binaries hav e been published ( Petrovic et al.|2005|[Cantiello| 
[et al. [[2007) 1. Our improvements to MESA now allow for the 
calculation of differentially-rotating binary stars. 

The rapid expansion of extra-solar planet research has led 
to a revival of interest in the detailed properties of stars probed 
through space-based brightness variability studies and radial 
velocity measurements. Stellar properties can be derived from 
measurements of the radial and non-radial oscillation modes 
of a star, but this requires the accurate and efficient compu¬ 
tations of mode frequencies and their eigenfunctions enabled 
by the coupling of GYRE with MESA. 

There are many ways M > 8 Mq stars can end their lives 
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(e.g,, Woosley et ar]|2002 Smart^2009| [Meynet et al.||2010 


Langer||2012| |Nomoto et al.||2013| |Smith||2014|l. Some be 

come electron capture supernovae; others collapse with most 
of their extended envelope intact and yield a Type II super¬ 
nova; others can lead to pair instability; and some have en¬ 
velopes thin enough to allow a jet to break through and appear 
as a long gamma-ray bur st ( MacFadyen & Woosley] T999[ 
Woosley & Bloom 2006 [Gehrels et al. 2009| l. There is a 


pressing need throughout the stellar community to routinely 
explore this entire mass range with new supernova progenitor 
and explosion models. The observational facilities discussed 
above have found explosions that indicate lar ge amounts o f 
mass are lost within a few years of explosion ( |Smith||T014 l; 
some show evidence o f optically thick w inds present at the 
moment of explosion ( jOfek et al.||2014|l , while others have 
yet to be securely identihed with a specific core collapse sce¬ 
nario. These mysteries, coupled with the community’s call for 
new yields from massive stars for galactic chemical evolution 
studies motivate the development of implicit shock hydrody¬ 
namics and explosions with large, in situ reactions networks 
in MESA. 

The paper is outlined as follows. Section describes the 
new capability of MESA to evolve binary systems. Section 
[^discusses the new non-adiabatic pulsation capabilities re¬ 
sulting from fully coupling to the GYRE software instrument. 
Section describes the improvements to accommodate im¬ 
plicit hydrodynamics with shocks. New capabilities for ad¬ 
vanced burning and X-ray bursts with large, in situ reaction 
networks are described in Section |5] In Section we model 
the pre-supemova evolution of massive stars and combine the 
implicit hydrodynamics module and the new capabilities for 
advanced burning to probe the nucleosynthesis and yields of 
core-collapse supernovae. Section [7] discusses the improve¬ 
ments for a more robust and efficient treatment of mass ac¬ 
cretion. Sectionl^presents a new option for an on-the-fly cal¬ 
culation of weak reaction rates and their application to the 
Urea process and accretion-induced-collapse models. Section 
1^ presents improvements in the physics implementation of 
particle diffusion by including radiative levitation and pushing 
diffusion methods into the strongly coupled, electron degen¬ 
erate regime. In Section we discuss improvements to the 
MESA software infrastructure, highlighting bit-for-bit consis¬ 
tency across operating systems and compilers. We conclude 
in Section[TT]by noting additional improvements to MESA are 
likely to occur in the near future. Important symbols are de- 
hned in Table [T] We denote components of MESA, such as 
modules and routines, in Courier font, e.g., evolve.star. 


Table 1 

Variable Index. 


Name 

Description First Appears 

a 

Orbital seperation 

2.1 

A 

Atomic mass number 



a 

Fine structure constant 



c 

Speed of light [z 


e 

Specific thermal energy 

4.4 

ri 

Wind mass loss coefiicient 

srr 

g 

Gravitational acceleration 

53 

G 

Gravitational constant 

2.1 

r 

Coulomb coupling parameter 


y 

/ 

Moment of inertia 

r 


K 

Opacity 



L 

Luminosity 

T 


A 

Reaction rate 



m 

Lagrangian mass coordinate 


zr 

M 

Stellar mass 

4 



Table 1 — Continued 


Name 

Description First Appears 

Ml 

Donor mass 


M 

M2 

Accretor mass 


>.A 


Mean molecular weight 


rr 

N 

Neutron number 


rr 

(Jj 

Dimensionless eigenfrequency 


).i 

n 

Angular frequency 


LA 

Q 

Nuclear rest mass energy difference 


rr 

p 

Pressure 



q 

Fractional mass coordinate 



qi 

Mass ratio, Mj /M 2 


23 

92 

Mass ratio, M 2 /Mi 


T3 

r 

Radial coordinate 

2 

1.1 

R 

Stellar radius 


P 

Baryon mass density 

2!: 


s 

Specific entropy 


i 

<J 

Oscillation eigenfrequency 



t 

Time 



T 

Temperature 


1.1 

T 

Timescale 


LA 

V 

Velocity 


rr 

X 

Baryon mass fraction 


h 

Y 

Molar abundance 


5 

z 

Gravitational redshift 

r 

53 

z 

Atomic number 



ffMCT 

Mixing length parameter 


13 

Cp 

Mass specific heat at constant pressure 



Xp 

(cJ \nP/dlnp)T 


7 

XT 

(a In P/a In T)p 


7 

dm 

Mass of cell 

2 

13 

dm 

mass associated with cell face 

2 

1.3 

dq 

Fractional mass of cell 



St 

Numerical timestep 

2 .; 

>3 

AM 

Change of stellar mass in one step 


/ 

Vad 

Adiabatic temperature gradient (<9 In r/<9 In P)ad 


7 

Vj 

Stellar temperature gradient d In T/d In P 


7 

Ef 

Fermi energy 

( 

d3 

^grav 

Gravitational heating rate 

2 

i.y 

fv 

Neutrino energy loss rate 

2 

i.y 

^nuc 

Nuclear energy generation rate 

2 

w 

^visc 

Viscous heating rate 

2 

1.4 

^visc 

Artificial dynamic viscosity coefficient 

2 

1.1 

/ov 

Convective overshoot parameter 


rr 

^visc 

Viscous acceleration 

2 

1.3 


Radiative acceleration 

y 

1.3 

Ti 

First adiabatic exponent {d In P/d\n p)ad 

233 

Hp 

Pressure scale height 


T3 

J 

Rate of change of angular momentum 


l.l 

dovh 

Orbital angulai' momentum 


LA 

kB 

Boltzmann constant 


1.1 

'^ion 

Mean inter-ion spacing 


i 

mp 

Proton mass 


13 

A/acc 

Accreted mass accumulated, Mt 


m 

Me 

Mass of unmodeled inert core 


:E 

^ign 

Mace at time of nova runaway 


73 

Mndd 

Eddington accretion rate 


lA 

Pn 

Electron chemical potential 


1.1 

^ion 

Ion number density 


i 

^osc 

Linear oscillation frequency 


l*! 

^’orb 

Orbital period 


1.1 

Gvisc 

Artificial viscosity energy 

2 

rz 

Rrl 

Roche lobe radius 


13 

Rd 

Debye radius 

y 

1.1 

Pc 

Central baryon mass density 


[6 

f 

GR con'ected time for observer at infinity 



Tc 

Central temperature 



TeS 

Effective temperature 

2 ., 

13 

^■dCC 

Timescale to accrete outer star layer 


£ 

rMLT 

Convective timescale 

2 

13 

"^OSC 

Oscillation e-folding time 


1.1 

Tsync 

Tidal synchronization timescale 


1.4 

Tth 

Thermal timescale of outer star layer 


£ 

rMLT 

Convection time scale 

2 

13 

V 

Time centered velocity 

2 

1.1 

Yc 

Electrons per baryon 
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2. BINARIES 

MESAbinary is a MESA module that uses MESAstar to 
evolve binary systems. It can be used to evolve a full stel¬ 
lar model plus a point mass companion or to simultaneously 
evolve the structure of two stars. It optionally allows the mod¬ 
eling of systems including stellar rotation, assuming the axis 
of rotation of each star to be perpendicular to the orbital plane, 
accounting for the effects of tidal interaction and spin-up 
through accretion. The impleme ntation of MESAbinary ben- 
efits from early contributions by |Madhusudhan et al.| (|2006]l 
and |Lin et al.] P011| l who modeled mass transfer frotn a star 
to a point mass. 

Here we provide an overview of the modelled physical 
processes for circular binary systems and describe the tests 
against which we validate MESAbinary. 


2.1. Initialization of a Circular Binary System 

A binary system is initialized by specifying the components 
and either the orbital period Porb or separation a. Each com¬ 
ponent can be a point mass or a stellar model. The initial 
model(s) are provided by a saved MESA model or a zero-age 
main-sequence (ZAMS) specification. Eor stellar models in¬ 
cluding rotation, the initial rotational velocities of each com¬ 
ponent can be explicitly defined, or set such that the star is 
synchronized to the orbit at the beginning of the evolution. 
The orbital angular momentum of the system is 

I Go 

7orb=MiM2J--— , (1) 

y Ml + M 2 

where Mi and M 2 are the stellar masses. Evolution of Mi, M 2 , 
and 7oib is used to update a using Equation Q. Masses can 
be modified both by Roche lobe overflow (RTOE) and winds. 
The total time derivatives of the component masses are given 
by 

Ml - ^l.w + A^rlOF, M2 - M2,w - fmtMRbOF, (2) 


of the Hulse-Taylor binary pulsar over three d ecades (Weis 
berg & Taylor|20()5]l and of the double pulsar (Kramer et al. 


2006|l have tested the predicted effect from gerieral relativity 


to a high precision. The angular momentum loss from gravi¬ 
tational waves is 


. 32 llnCV^^ (MiM2f 

(Mi+M2)2/3- 

2.2.2. Mass Loss 

We assume the mass lost in a stellar wind has the specific 
orbital angular momentum of its star. Eor the case o f ineffi- 
cient mass transfer, angular momentum loss follows |Sober-| 
[man et 3T] ( |1997[ ), where fixed fractions of the transferred 
mass are lost either as a fast isotropic wind from each star 
or a circumbinary toroid with a given radius: 

Tml = [(47l,w + Q'mt47RLOF)47j + (M2,w + ^mtA^fRLOF)^^^ j 

a^ 2 ti 

^ {Ml + M2)2 T^b 

+ Tmtdmt-^RLOF sjGiMi + M2)a, (5) 

where amt, Pmu and (5mt are respectively the fractions of mass 
transferred that is lost from the vicinity of the donor, accretor 
and circumbinary toroid, and is the radius of the toroid. 
Ignoring winds, the efficiency of mass transfer is then given 
by /mt = 1 “ Q'mt “ Pmi “ ^mt- When accretion is limited to 
MRdd, efficiency of accretion is given by 

/mt ~ min(l — amt ~ (Sml ~ dmt, |4/Edd/4/RLOFl), (6) 

and the additional mass being lost is added to the PmiMRbOF 
term in Equation Q, i.e., it is assumed to leave the system 
carrying the specific orbital angular momentum of the accre¬ 
tor. 


where Mi is the donor mass and M 2 the accretor mass. The 
stellar wind mass loss rates are and M 2 ,w (see Paper I and 
Paper II) and Mrlof is the mass transfer rate from RLOE, all 
defined as negative. The factor /mt represents the efficiency of 
accretion and can be used to limit accretion to the Eddington 
rate Medd- 

2.2. Evolution of Orbital Angular Momentum 

To compute the rate of change of orbital angular momen¬ 
tum, we consider the contribution of gravitational waves, 
mass loss, magnetic braking, and spin orbit (LS) coupling 

Torb “ 7gr -F Tml "t" 7mb + /s , (3) 

from which the change in orbital angular momentum in one 
step is calculated as ATorb = /orbdf, where 6t is the timestep. 
Unless models with stellar rotation are being used, the /s term 
is equal to zero, and the contribution of the individual spins 
of each star is not directly considered. On the other hand, the 
/mb term implicitly assumes a strong tide that keeps the orbit 
synchronized. The simultaneous us age of /mb with stellar ro¬ 
tation is not consistent (see Section [2. 2. 4| i. We now describe 
how these terms are computed. 

2.2.1. Gravitational Wave Radiation 

Very compact binaries can experience significant orbital de¬ 
cay due to the emission of gravitational waves. Observations 


2.2.3. Spin Orbit Coupling 

Tidal interaction and mass transfer can significantly mod¬ 
ify the spin angular momentum of the stars in a binary sys¬ 
tem, acting as both sources and sinks for orbital angular mo¬ 
mentum. The impact spin-orbit interactions have on orbital 
evolution depends on the orbital separation and the mass ra¬ 
tio, with the effect being greater for tighter orbits and uneven 
masses. The corresponding contribution to /orb is computed 
by demanding conservation of the total angular momentum, 
accounting for losses due to the other /orb mechanisms and 
loss of stellar angular momentum due to winds. 

In a fully conservative system, the change in orbital angular 
momentum in one timestep is 6Jorb - SS i - 6S 2 , where 
6S 1 and 6S 2 are the changes in spin angular momenta. This 
needs to be corrected if mass loss is included, as winds take 
away angular momentum from the system. If S i.iost and S 2 ,lost 
are the amounts of spin angular momentum removed in a step 
from each star due to mass loss (including winds and RLOE), 


J\s - 



„ 4 / 1 ,w 

^ l,lost 

Ml 



(7) 


where the additional factor for the donor accounts for mass 
lost from the system, ignoring mass loss due to mass transfer. 
In the absence of RLOE this equation becomes symmetric be¬ 
tween both stars, as then = L 
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The form of Equation (j^ is independent of how tides and 
angular momentum accretion work, as it is merely a statement 
on angular momentum conservation. The details of how we 
model these processes and their impact on the spin of each 


component are described in Section 2.4 


2.2.4. Magnetic Braking 

The rotational velociti es of low mass star s are strongly cor¬ 
related with their ages (jSkumanich|| 1 972)1. This spin-down 
arises from the coupling of the stellar wind to a magnetic field. 
If the star is in a binary system and tidally coupled to the orbit, 
magnetic braking ca n provide a very efficient sink for orbita l 
angular momentum ( |Mestel |1968[ [Verbunt & Zwaan] 
lis effect following Rapp 


1981 


We implement this effect following Rappaport et al. 
who assumed the star being braked is tidally synchronized: 


=-6.82 X 10^^ ^1 [dyncm], (8) 


where in the simplest approximation ymb=4 ( [Verbunt 
|Zwaan||1981| l. A similar contribution from the accretor can 
be included. As tidal synchronization is assumed, this formu¬ 
lation is incompatible with the use of LS coupling. 

It is normally assumed that once a star becomes fully con¬ 
vective, the dynamo process that regenerates the field will stop 
working or at least behave in a qualitatively different manner. 
Similarly, magnetic fields in stars with radiative envelopes are 
of a significantly different nature than those seen in stars with 
convective envelopes, and there is no simple way to predic t 
even the presence of magnetism (jDonati & Landstreet|2009|l. 
By default MESAbinary only accounts for magnetic braking 
as long as the star being braked has a convective envelope and 
a radiative core, though the process might still operate outside 
of these conditions. 


2.3. Mass Transfer from RLOF 

Close binary stars are defined as systems tight enough to 
interact through mass transfer, with the most important mech¬ 
anism being RLOF. This process is commonly modeled in 
ID by considering t he spherical-equ ivalent Roche lobe radius 
Rrl of each object (|Eggleton|1983|l 


R 


QA9q 


2/3 


0.6q^/^ + 


ln(l 




a, 


(9) 


where j is the index identifying each star, qi - Mi I M 2 and 
q 2 = M 2 /Ml. This fit is correct up to a few percent for the full 
range of mass ratios, 0 < qj < oa. Mass transfer occurs then 
when the radius of a star approaches or exceeds Rrl. Depend¬ 
ing on several factors, once a star begins RLOF the ensuing 
mass transfer phase can proceed on a nuclear, thermal, or dy¬ 
namical timescale. 

The stability of mass transfer is normally understood in 


terms of mass-radius relationships (e.g.. Tout et al. 1997 
jSoberman et al.|1997| . 


Here, gives the radial response of the donor to mass loss 
when it happens slowly enough for the star to remain in ther¬ 
mal equilibrium. When mass loss proceeds on a timescale 
much shorter than the thermal timescale of the star, but still 
slow enough for the star to retain hydrostatic equilibrium then 
the radial response will be given by ^ad- The dependency of 
the Roche lobe radius on mass transfer is encoded in ^rr. In 
general = dlnRi/dlnMi is a function of Mrlof, so re¬ 
quiring C, - ^RL will determine the value of Mrlof- If an 
overflowing star satisfies feq > ^rl, then it can remain inside 
its Roche lobe by transferring mass while retaining thermal 
equilibrium. If on the contrary ^ad > ^rl > mass trans¬ 
fer will proceed on a thermal timescale, while for the extreme 
case /^RL > /fad the star will depart from hydrostatic equilib¬ 
rium and the process will be dynamical. MESA cannot model 
common envelope or contact binaries. 

MESAbinary provides both explicit and implicit methods to 
compute mass transfer rates. An explicit computation sets the 
value of Mrlof at the start of a step, while an implicit one 
begins with a guess for Mrlof and iterates until the required 
tolerance is reached. The composition of accreted material is 
set to that of the donor surface, and the specific entropy of 
accreted material is the same as the surface of the accretor. In 
models with rotation the specific ang ular momentum of ac¬ 
creted material is described in Section l24l 

2.3.1. Explicit Methods 

MESAb inary implem ents two mass transfer schemes: the 
mod el of|Ritter|(|1988]l whic h we refer to as the Ritter scheme 
and [Kolb & Ritterj ( j 19^0) 1 which we refer to as the Kolb 
scheme. We use the mass ratio q 2 consistent with the Ritter 
scheme. 


Ritter scheme: — Stars have extended atmospheres therefore 
RLOF can take place through the LI point even when Ri < 
Rr 


'RL,1- 


Ritter (1988 1 estimated the mass transfer rate for this 


case as 


i R\ — Rrl 1 \ 

Mrlof = -Mq exp ‘ , (13) 

\Hp,ilyiq2)l 

where Hp i is the pressure scale height at the photosphere of 
the donor and 


Mo 


2n 

exp(l/2) 


Ffqi) 


4l.i 


GMi 


nip^ph) 


(14) 


where mp is the proton mass, Tetr is the effective temperature 
of the donor, and jUph and pph are the mean molecular weight 
and density at its photosphere. The two fitting functions are 


( 0 ^ 2 ) = 1-23-F 0.5 log g' 2 , 0.5<q'2^10. (15) 


and 


- j 

1 dlnRi \ 

(10) 

yiqi) = 

[dlnMi 

1 d\nRi \ 

(0.954 +0.025logg 2 - 0.038(log^2)^ 0.04 < ^2 < 1 
\0.954 + 0.039 log ^2 + 0.114(log^2)^ 1 < ^2 < 20. ^ ^ 




(11) 

[dlnMi 

Outside the ranges of validity Fi{q 2 ) and y{q 2 ) are evaluated 


/"rl = 

£/ln/?RL,i 

(12) 

using the value of q 2 at the edge of their respective ranges. 

InMi 
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Kolb scheme: — [Kolb & Ritter] ( |1990| ) extended the Ritter 
scheme in order to cover the case R\ > Rrl,i according to 


A^^rlof - -Mq - 2nF\ {q 2 ) 


Ri 


GMi 


: r Tf 


r, +1 


(ri+i)/(2ri-2) 


\mpfi 


1/2 


dP (17) 


where Ti is the first adiabatic exponent, and Pph and Prl are 
respectively the pressures at the photosphere and at the radius 
for which ri = Rrl.i- 

2.3.2. Implicit Methods 

Explicit schemes exhibit large jumps in Mrlof unless the 
timestep is severely restricted. Therefore, if one needs accu¬ 
rate values of Mrlof and stellar radius, this requires use of an 
implicit scheme. Implicit schemes also allow the calculation 
these quantities when there is no general closed form formula 
for Mrlof- 

These implicit methods use a bisection-based root solve to 
satisfy |/(Mrlof)I < ^ at the end of the step, where ^ is a 
given tolerance. The implicit schemes are then defined by the 
choice of the function /(Mrlof)- For the Ritter and the Kolb 
scheme the function is chosen as 


/(.Mrlof) 


Mend - Mrlof 


(18) 


with Mend being the mass transfer rate computed at the end of 
each iteration. 

A different implicit method is also provided. In this case, 
whenever the donor star overflows its Roche lobe the implicit 
solver will adjust the mass trans fer rate until Ri - Rrl.i 
within some tolerance (see e.g., Whyte & Eggleton 1980[ 
Rappaport et al.|1982 1983|l. In this case 


/(Mrlof) 


2(Ri - Rrl.i) 
f^RL.l 




(19) 


and if Mrlof is below a certain threshold and /(Mrlof) < 
then the system is assumed to detach and Mrlof is set to zero. 

2.4. Effect of Tides and Accretion on Stellar Spin 


To model tidal interaction we adjusted the model of Hut 
(|1981|l to include the case of differentially rotating stars. The 
time evolution of the angular frequency for each component 


dTlc i ^orb 


dt 


‘ sync,/ 


* sync,/ 


(-) (20) 


i^jrgJP \TI,j\ a 


where j = 1,2 is 
gular frequency at 


the index of each star, Q, y is the an- 
the face of cell i towards the surface, 
rG = Ijl(MjRj) is the radius of gyration (with Ij being the 
moment of inertia of each star), and the ratio of the apsidal 
motion constant t o the viscous dissipat ion timescale, (A:/r)c,; 


is computed as in Hurley et al. 


( 2002| l. Similarly to lDetmer s 
ant Tsync,; and Qorb through 


et al. (20081, we assume constant Tsync,; and iiorb ttirougti a 
step and therefore = [1 - exp(-hf/Tsync,;)](f^oib “ 

This extension of Hut’s work to differentially rotating stars 
is not formally derived but merely applies his result for solid 
body rotato rs independently to each shell. The formulation of 
Hut (19811 can be recovered from Equation ([20|), by forcing 


solid body rotation with a large diffusion coefficient for an¬ 
gular momentum throughout the star. In reality tides would 
act mostly on the outer layers, and whether the core synchro¬ 
nizes or not depends on the coupling between the core and the 
envelope. 

To compute the specific angular momentum carried by ac¬ 
creted material, we consider the possibi lity of both ballistic 
and Keplerian disk mass transfer (e.g., [Marsh et al.||2004 
de Mink et al.|[2013|l. To distinguish which occurs, we 


compar e the minimurn distance of approach of the accretion 
stream ( Lubow & Shu||l975j [Ulrich & Burger 1976p 


f^min = 0.0425a (g'2 + ql) 


1/4 


0.0667 < ^2 < 15 (21) 


to the radius of the accreting star. When outside the range of 
validity, R^in is computed using the value of q 2 at the respec¬ 
tive edge. Accretion is assumed to be ballistic whenever R 2 > 
Rmin and the specific angular momentum is (1.7GM2Rrnin)*^^- 
When R 2 < Rmin the specific angular momentum is taken as 
that of a Keplerian orbit at the surface (GM 2 R 2 )'^^- 

2.5. Treatment of Thermohaline Mixing in Accreting Models 

In stars with radiative envelopes accreted material with a 
high mean molecular weight is expected to mix inwards due 
to thermohaline mixin g, a process that is very sensitive to 
the ju-gradien t (see e.g., Kippenhahn et al.|1980[ Cantiello & 
|Langer|2010| l. Thetmohaline mixing is included in MESA (see 
Paper 1). However, as mass with homogeneous composition 
is added during the accretion process, a jump is produced at 
the boundary between new and old material. MESAstar com¬ 
putes mixing coefficients explicitly at the start of each step, 
so this results in thermohaline mixing only operating near this 
boundary, leading to unphysical compositional staircases. To 
avoid this issue, we artificially soften the composition gra¬ 
dient in the outer (A(7)iarge fraction of the star by mass. We 
do this starting at the surface and homogeneously mixing in¬ 
wards a region of size (Aq’jsmaii. Then, moving towards the 
center, the process is repeated at each cell while linearly (with 
respect to mass) reducing the size of the small mixed region 
such that it is zero after going (A^)iarge inwards. All the binary 
models where the accretor is not a point mass are calculated 
using (A^)iarge = 0.05 and (A^)smaii = 0.03. 

2.6. Numerical Tests 

Here we describe tests designed to va lidat e the implementa¬ 
tion of the physics described in Section 12.21 We check orbital 
evolution in the presence of gravitational waves and mass loss 
by comparing to analytical solutions. We also verify total an¬ 
gular momentum conservation in calculations that include the 
physics of tides and spinup by accretion. To test for the ther¬ 
mal response of stellar models undergoing mass transfer, we 
compare MESAbinary results to those from the STARS code 


( Eggleton] 1971 Pols et al. 1995t Stancliffe & Eldridge|2009 1. 


2.6.1. Gravitational Wave Radiation 

If gravitational waves are the only source of angular mo¬ 
mentum loss and the masses of each component remain con¬ 
stant, Equation © can be i ntegrated to o btain the time evo¬ 
lution of orbital separation (|Peters||1964|l. We model a sys¬ 
tem consisting of a 0.5 M© star and a 0.8 M© point mass with 


^ ^ Note that there is a small typo in the fit given by | U lrich & Burger|(l976 . 
The corrected fit given here fits the results of lLubow^^EuHl^VS^ to the 4% 
accuracy claimed by |Ulrich & Burg^jl976| . 
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log (Toff/K) 

Figure 1. Evolution in the HR diagram for a 2.5 M© star transfening mass to 
a 1.4 M© point mass, assuming a mass transfer efficiency of 1 %. Symbols are 
shown at zero-age main sequence (ZAMS) and terminal-age main sequence 
(TAMS), together with parts of the track where RLOF is occuning. The inset 
shows evolution from ZAMS up to the beginning of the first phase of mass 
transfer. 



A/i [Me] 

Figure 2. Evolution of mass transfer rate from a 2.5 M© to a 1.4M© point 
mass, assuming a mass transfer efficiency of 1%. The upper panel shows the 
difference between the computed orbital separation and the analytical solu¬ 
tion while the bottom one displays the evolution of the mass transfer rate, 
using two different schemes. 


a - 2 Rq. We ignore all effects on the evolution of orbital an¬ 
gular momentum except its loss due to gravitational waves. 
In 3.5 Gyr the orbital separation of this system reduces to 
a = 1.3 Rq, at which point the 0.5 Mq star begins mass trans¬ 
fer. We terminate the run at the onset of RLOF. The maximum 
error in a is 0.35% relative to the analytical result. 


2.6.2. Inejficient Mass Transfer 

An analytical expression for the evolution of orbital sepa¬ 
ration can be derived if inefficient mass transfer i s the only 
contribution to the angular momentum evolution (Tauris & 


van den Heuvel 


20061. We model a 2.5 Mq main sequence 
star together with a 1.4 Mq point mass with an initial orbital 
separation of 10 Rq. We choose a^t = 0.03, jSmt = 0.95, 
dmt = 0.01 and = 2, which give a low mass transfer effi¬ 
ciency of /jnt = 0.01. Such a system is representative of the 
evolution of an intermediate mass X-ray binary (IMXB). The 
model initiates mass transfer just after the end of the main 
sequence, interrupting the evolution of the star through the 
Hertzsprung gap and producing a low mass white dwarf (WD) 
(Mhb = 0.289 Mq) with a small amount of hydrogen on its 
surface. As the WD evolves to the cooling track, it experi¬ 
ences several hydrogen flashes, one of them strong enough to 
produce an additional phase of RLOF (see Figure [^. 

Figurej^shows that MESAbinary computes the orbital evo¬ 
lution to a precision of a few parts in 10^. We rain this sys¬ 
tem using both the Ritter and the Kolb implicit schemes to 
display that under some circumstances the precise choice of 
mass transfer scheme does not play a big role in the evolution. 


2.6.3. Spin Orbit Coupling 

We now test angular momentum conservation by ignoring 
all the mechanisms that remove angular momentum from the 


binary system. For this purpose we model an 8 Mq + 6 Mq 
binary with rotating components and an initial orbital period 
of 1.5 days. Due to the short orbital separation we assume 
the initial spin periods of the two stars are equal to the or¬ 
bital period. The primary undergoes RLOF during the main 
sequence, initiating a phase of mass transfer on a thermal 
timescale. After transferring just 0.3 Mq the accretor also hlls 
its Roche lobe, producing a contact system. At this point we 
terminate the evolution. 

Figure]^ shows that spin angular momentum in both com¬ 
ponents increases during the pre-interaction phase, which is 
due to both stars expanding on the main sequence while re¬ 
maining tidally locked. During Roche lobe overflow, the sec¬ 
ondary is rapidly spun-up, reaching nearly 80% of critical ro¬ 
tation before contact. The calculation of total angular momen¬ 
tum requires the summation of different contributions (orbital 
angular momentum and spin of both components). Therefore 
the maximum accuracy to which we can conserve angular mo¬ 
mentum is limited by rounding errors. Figure shows that 
conservation of angular momentum in the run is very close to 
machine precision. 


2.6.4. Thermal Response to Mass Loss 

The fate of binary systems depends largely on the precise 
value of M during an interaction phase, which depends on the 
thermal response of the donor star to mass loss. For WDs 
there is a limited range of accretion rates for stable hydrogen 
burning (Nomoto et al. 2007 Shen & Bildsten|2007| l. In main 
sequence binaries the evolution of the accretor radius depends 
on the mass transfer rate, and expansion durin g the interaction 
phase can lead to contact or even a merger (IWellstein et al.| 

[200T] l. 

We calculated an 8 Mq -h 6.5 Mq binary system with an ini¬ 
tial orbital period of 1.5 days using both MESAbinary and 
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Figure 3. Angular momentum evolution in an 8 M© + 6 M© binary with 
an initial orbital period of 1.5 days. Left panels show the evolution before 
the onset of RLOF, while right panels display evolution from the beginning 
of RLOF until contact, when both components fill their Roche lobe. The 
fractional error in the total angular momentum is plotted in the bottom panel 
and is of order machine-precision. 


STARS. To minimize the modeling differences and focus on 
the thermal response of both components, we use an ex¬ 
tremely simplified model that ignores internal mixing (includ¬ 
ing convective mixing). Under these conditions, the more 
massive star quickly depletes its central hydrogen and be¬ 
gins shell hydrogen burning, reaching RLOF and undergoing 
a phase of mass transfer on the thermal timescale. The result¬ 
ing mass transfer rates are shown in Figure The agree¬ 
ment is very good, despite mass transfer rates being com¬ 
puted in slightly different ways. Masses at detachment show 
a small difference, with the MESAbinary model terminating 
mass transfer when Mi = 0.952 Mq while the STARS cal¬ 
culation when M] = 0.935 Mq. The hgure also shows the 
change in radius of the accreting star, with two prominent 
peaks at R 2 /R 0 = 4.84,5.34 for MESAbinary and 4.82,5.28 
for STARS. The larger radius of the MESAbinary model is 
likely associated to the slightly higher mass transfer rates. 


2.7. Period Gap of Cataclysmic Variables 

Although cataclysmic variables (CVs) span a wide range 
of periods, observations show a lack of systems in the rang e 
2 h < Porb < 3 h (see, for instance, [Gansicke et al.|[2b09| l. 
Such a feature is commonly explained by having an angular 
momentum loss mechanism “turn off” or become inefficient 
at some point. The m ost popular model for s uch a mechanism 
is magnetic braking ( [Rappaport et al.|1983[ ), as the magnetic 
held of the donor is aissumed to change quickly when the star 

loses enough mass to become fully convective. _ 

We now compare to the results of|Howell et al.|(|2001|l, who 
performed a population synthesis study to explore in detail the 
standard scenario involving magnetic braking. In Figure|^we 
show the evolution of mass transfer rates and orbital periods 
for a set of CV models with different component masses and 



Figure 4. Mass transfer rate and accretor radius as computed by MESA and 
STARS for an 8 M© + 6.5 M© binary with an initial orbital period of 3 days. 
All internal mixing processes (including convective mixing) ai‘e turned off in 
the calculations. 



Figure 5. Evolution of CV models under the effect of magnetic braking and 
gravitational wave radiation. For each track the label gives the donor mass, 
the WD mass, and the initial orbital period respectively. The grey ban d shows 
the observed period gap for CVs. These results reproduce figure 1 in|Howell| 

FnrjpooT} - - 


orbital periods. We run all models using - 1 and ymb = 3 
and magnetic braking is turned off when the donor star be¬ 
comes fully convective. As an example the system with a 
0.9 Mq donor (left panel in Figure]^ experiences a first phase 
of mass transfer induced by magnetic braking between 10^ ' 
and 10^'^ yr, a non-interacting phase (the gap) between 10* ^ 
and 10^ ^ yr, and a subsequent phase of mass transfer dom¬ 
inated by gravitational wave radiation, reaching a minimum 
orbital period of ab out 1 hour at 10^-^ yr . As a comparison, 
for the same model |Howell et al.| (|2001|l obtain a first phase 
of mass transfer between 10^ ^ and 10® “* yr, the gap occurs be¬ 
tween 10® “* and 10® ® yr and a period minimum is reached at 
10® “* yr. Figure^ shows that our CV models spend most time 
away from the ^served period gap. 


2.8. Evolution of Massive Binaries 
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Figure 6 . Evolution of a 16 Me + 14Mo system with a 3 d ay initial orbital 

S ieriod. MESAbinary models are compared to the results of|Wellstein et al.| 
200 1(, which were calculated using the STERN code. The terms primary and 
secondary are used throughout the evolution to describe the initially more 
massive and the less massive components, respectively. For each component 
in the MESAbinary model, squares mark the ZAMS and the depletion of the 
indicated nuclear fuel in the core. 


In massive stars, binary interactions have dramatic effects 


on the evolution of both components. Kippenhahn & Weigert 
( 1967[) introduced the term “case A” to refer to a mass trans 


fer phase occurring in systems tight enough such that RLOF 
starts during the main sequence. This results in a large amount 
of mass being transferred on a thermal timescale, followed by 
a phase of mass transfer that proceeds on the nuclear timescale 
until the end of core H-burning. An additional phase of 
thermal-timescale mass transfer then follows (the so-called 
“case AB”), which strips the donor and produces an almost- 
naked helium star. 

Here we show that MESAbinary can calculate the evolution 
of massi ve interacting binaries . We reproduce one of the mod¬ 
els from|Wellstein et al.|(|20011l, a 16 Mq + 14Mo system with 
an initial period of 3 days, using the same semiconvection effi¬ 
ciency of Q-sc = 0.01. As shown in Figures|^and|^this system 
experiences case A and AB mass transfer, and the accretor be¬ 
comes a blue supergiant after core hydrogen depletion. The 
accretor depletes carbon before its donor. 

Figure 1^ illustrates the prevalence of both thermohaline 
mixing am semiconvection in the accreting star. Newly ac¬ 
creted material is efficiently mixed inwards by thermohaline 
mixing. On the other hand the p-gradient formed before in¬ 
teraction prevents the convective core from growing, with the 
efficiency of semiconvection controlling whether or not the 
star rejuvenates. Due to the choice of inefficient semicon¬ 
vection, the core remains small, preventing the star from be¬ 
coming a red supergiant. The star accretes a large amount of 
CNO-processed and helium-rich material. After being mixed 
through the envelope this material results in the surface being 
nitrogen rich and carbon depleted, with a slight enhancement 
in helium. 


Figure 7. Kippenhahn diagram for the evolution of a 16Mo + 14 Mq system 
with a 3 day initial orbital period. Most of the pre-interaction phase is not 
shown in this figure. The upper plot shows the evolution of the donor, while 
the lower plot displays that of the accretor. 


2.9. Rotating Binaries and the Efficiency of Mass Transfer 

The efficiency of mass transfer plays a key role in close 
binary systems, but the processes by which material is lost 
from the system are not well-understood. In particular, when¬ 
ever an accreting star approaches Q/Qdit = 1, it is uncertain 
whether accretion can continue, one option being the devel- 
opment of a strong wind that prevents accretion (e.g., |Petro-| 
vie et al.|2005[ Cantiello et aH2007) l. Whenever Q/Qcrit ap- 
proaches one, we use an implicit method to iteratively reduce 
M 2 until this ratio falls below a threshold. 

Tides counteract the effect of spin-up from accretion. 
Whether or not an accreting object reaches critical rotation 
depends on the efficiency of tidal coupling. Here we model 
a I 6 M 0 -I- 15Mo binary system including differential stellar 
rotation, with an initial orbital pe riod of 3 days and a ssum- 
ing initial orbital synchronization. Danger et al.| ( |2003] l argue 
that turbulent processes in the radiative envelope can signif¬ 
icantly enhance tidal strength. They model the same system 
using the simple estimate for the synchronizat ion timescal e 
for a star with a convective envelope given by Zahn (19771, 
Tsync.f = 1 yr X qj{alRf. For our implicit modeling of stellar 
winds we use a threshold of (Q/Qcrit) = 0.99. 

Figure s hows that MESAbinary m odels using both the 
Zahnj (|1977| and |Hurley et al.| (|2002|l timescales for tidal 
coupling. These rnodels experience highly non-conservative 


pling. These rnodels experience highly 
phases of mass transfer, corresponding to the accreting star 
evolving very close to critical rotation. In particular during 
case AB mass transfer the accretor needs to switch from mass 
accretion to mass loss in order to remain sub-cr itical. As ex- 
pected, the system with the tidal timescale from |Zahn| (| 1 977)1 
has a significantly higher mass transfer efficiency, and during 
the first phase of RLOF it only experiences a brief period in 
which the accretor reaches c ritical rotation. This is in broad 
agreement with the model by|Langer et al.|p003|l. 
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Figures. Efficiency of mass transfer in a 16Mo + 15 Mq binary system 
in cluding differential r otation. The system is modeled with tides as described 
by |Hurley et al.|(2002) for radiative envelopes, and also with the simple tidal 
tiniescale given by |Ii^ahn|jl977). The upper panel shows the efficiency of 
mass transfer, the rniddle panel the angular frequency of each star in terms 
of its critical value, while the lower panel shows the evolution of M for both 
components. 


2.10. Description of a Binary Run 

MESAbinary performs each evolution step by indepen¬ 
dently solving the stracture of each component and the orbital 
parameters, using the same timestep 6t for each. This ap¬ 
proach differs from STARS, which simultaneously solves for 
the structure of both stars and the orbit in a single Newton- 
Raphson solver. Our choice to solve for each star separately 
gives a significant amount of flexibility and simplicity, as the 
examples in this paper demonstrate. 

The top-level algorithm for evolving a star is described in 
appendix B1 of Paper II. We modified this algorithm to sup¬ 
port the new implementation of binary interactions, which is 
described in detail in the MESA documentation. Additional 
timestep limits are imposed in MESAbinary that consider rel¬ 
ative changes between the radius and Roche lobe radius of 
both components, the total orbital angular momentum, the or¬ 
bital separation, and the envelope mass in the donor. 


3. PULSATIONS 

The study of stellar pulsations (also termed oscillations) 
offers unique insights into the interiors of stars ( |Aerts et al.| 
2010|l. In some classes of star (e.g., solar-type, red giant), the 


stochastic excitation of hundreds of oscillation modes, typi 
cally by convective motions, allows remarkably detailed mea¬ 
surements tobemudeoftheinteri or, including nuclear bum 


ing state (Bedding et al. 20111 and internal rotation (Beck 


et al.||2012|l. In other classes (e.g., classical Cepheid, f 
Cephei, 6 Scuti and y Doradus pulsators), modes are instead 
excited by linear instabilities, most often linked to opacity 
variations in the envelope (the k mechanism). In these lat¬ 
ter objects, typically too few modes are excited for detailed 
asteroseismic analysis to be feasible; nevertheless, mapping 
out the regions of the theoretical HR diagram where the in¬ 
stabilities are expected to operate, and then comparing these 
instability strips against observational surveys, can often lead 
to new science. 

Paper II introduced the astero extension to MESAstar, 
which permits on-the-fly refinement of stellar model param¬ 
eters in order to fit a set of observed oscillation frequencies 
and spectroscopic constraints. Subsequent improvements to 
the a stero capabili t ies inc lude frequency correction recipes 
from [Ball & Gizon] (|20141l; im plementation of the downhill 
simplex (|Nelder & Mead|l%5|l and NEWYUO (|Powell|2004]l 
algorithrnsTor^plnmimization; parameter optimization using 
only spectroscopic constraints (e.g., Teff and surface gravity); 
and coupling to the GYRE oscillation code, as an alternative 
to the ADIPLS code ( Christensen-Dalsgaard|20d8l ) used in the 
original implementation. 

GYRE calculates the normal-mode eigenfrequencies cr of a 
stellar model by solving the system of linearized equations 
and boundary conditions governing small periodic perturba¬ 
tions (cx exp[icrf]) to the equilibrium state. It is based on a 
novel Magnus Multiple Shooting (MMS) numerical scheme 
which is robust and accurate, and makes full use of all avail¬ 
able processors on multicore computer architectures. The 
MMS scheme and the initial release of the code, which fo¬ 


cuses on adiabatic pulsations, is described in Townsend & 
Teitler|(|2013)l; extensions to the code to support non-adiabatic 
pulsations are described in Goldstein & Townsend (2015 i. 

MESAstar couples to GYRE via two mechanisms^ Loose 
coupling is achieved simply by MESAstar writing models 
out to disk, and GYRE subsequently reading these models in; 
we use this process below to map out massive-star instability 
strips. Tight coupling removes the intermediate disk usage, 
by handling all communication between MESAstar and GYRE 
in-memory; this permits fully closed-loop calculations, where 
the changes in the oscillation eigenfrequencies of an evolving 
stellar model are used to guide the further evolution of the 
model. Tight coupling allows GYRE to function as an alter¬ 
native to ADIPLS in the astero extension, and opens up the 
possibility of other kinds of novel calculations, such as the 
automated location of instability-strip boundaries. 

3.1. Massive-Star Instability Strips 

As an illustration of a large-scale calculation using 
MESAstar and GYRE loosely coupled. Figure 1^ plots the in¬ 
stability strips for massive stars on and near the upper main 
sequence, for oscillation modes with harmonic degrees £ - 
0-3. These strips are based on a set of 182 evolutionary 
tracks, each extending from the ZAMS across to a red limit at 
log(7’eff/K) = 3.75, with 101 tracks spanning the initial mass 
range 2.5 Mq < M < 25 Mq in uniform logarithmic incre¬ 
ments, and the remaining 81 tracks spanning the mass range 
6Mo < M < 10 Mq in uniform linear increments (the latter 
set is designed to adequately resolve the “fingers” discussed 
below). OPAL o pacity tables are used with the proto-solar 
abundances from |Asplund et al.| (|20091l, and for simplicity 
we neglect any rotation or mass loss. Convection is modeled 
with a mixing-length parameter (Tmlt = L5 and an exponen¬ 
tial overshoot parameter /ov = 0.024, and the Schwarzschild 
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Figure 9. Instability strips for i = 0-3 oscillation modes in the upper pai't of the Hertzsprung-Russell diagram. Separate strips are shown for the y? Cephei (tj > 1) 
and slowly pulsating B-type (SPB; tu < 1) classes of pulsating stars. The ZAMS and red edge of the main sequence (REMS) are shown for reference, as are 
evolutionary tracks for models with selected masses (labeled in solar units along the ZAMS). The red edges of the post-MS SPB strips are drawn with a dotted 
line, indicating that the positioning of these edges is an artifact of our numerical procedure. 


stability criterion is assumed. 

We select points i - h, ■ ■ ■ along each of the 182 tracks 
(where i is the timestep index; see section 6.4 of Paper I), 
chosen so that corresponds to the ZAMS, 


Z 

;'=/i 


friog(7’eff,,+ l/7’efr,,) 

2 

log(L,+i/L,) 

n 

[i At- 



I 


I (22) 


across the (/i.ia) pair, and similarly for subsequent pairs. 
Here, At and At are dimensionless weights which control 
the spacing of points in effective temperature and luminos¬ 


ity; we adopt the values 0.004 and O.Oll, respectively, for 
these weights. At the selected points, GYRE searches for un¬ 
stable oscillation modes with the harmonic degrees consid¬ 
ered. First, GYRE solves the adiabatic oscillation equations to 
find eigenfrequencies falling in the range extending from 
the asymptotic frequency of the gravity (g) mode with radial 
order n - 400, up to the asymptotic frequency of the pressure 
(p) mode with radial order « = lO. Each (Tad is then used as 
an initial guess in finding a corresponding eigenfrequency cr 
of the full non-adiabatic oscillation equations. The real and 
imaginary parts of cr give the linear frequency Vosc and the 
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log(Teff/K) 


Figure 10. The I = \ dimensionless frequency spectrum of an 8.5 M© stellar 
model as it evolves from the ZAMS to the REMS. Blue (orange) dots indi¬ 
cate which modes are stable (unstable); selected modes are labeled along the 
left/bottom edge using their classification. 



log(Teff/K) 


Figure 11. Instability strips for dipole (/’=!) oscillation modes in the upper 
part of the HR diagram, hut calculated using OP rather than OPAL opacities 
(cf. Figure]^. 


growth e-folding time Tqsc of a mode: 


_ Re(o-) _ _ 1 

2ji ’ Im(cr) 


(23) 


If Tosc is negative, the mode is damped. 

Separate strips are shown in Figure [^for regions exhibiting 


unstable modes with Re(aj) > 1 and Re(cu) < 1, where 


0 ) - cr 


R3 


(24) 


is the dimensionless eigenfrequency; these correspond, re¬ 
spectively, to the y6 Cephei and slowly pulsating B-type (SPB) 
classes of pulsating stars. In (3 Cephei stars during the MS 
phase, p and g modes with periods of a few hours and ra¬ 
dial orders n 1-3 are excited by a ic mechanism operat¬ 
ing on the iron opa city bump situated in the outer envelope 
at log(7’/K) ai 5.3 dCox et al.||1992 [Dziembowski & Pami-| 
atnykh|[T993| l. In SPB stars during the MS phase, g modes 
with periods of a few days and r adial orders n 20-50 ar e 
excited by the same mechanism ( [Dziembowski et ^|1993 l. 
For masses M > 9 Mq the strips for both classes of stars ex¬ 
tend into the post-MS domain. During this phase, unstable 
modes couple with g modes trapped near the boundary of the 
inert helium core. In the case of the SPB stars this leads to 
very high overall radial orders, n > 100, and ultimately limits 
our ability to follow the instability strips all the way to the red 
edge (our calculations are restricted to n < 400 for computa¬ 
tional efficiency reasons). Hence, in Figurel^we plot the red 
edges of the post-MS SPB strips with dotted lines, to highlight 
that these are not the true red edges. 

Allowing for differences in adopted abundances and other 
modeling parameters, the instability strips plotted in Figure]^ 
are in ge neral agreement with those published in the litera¬ 
ture (e.g.,|Pamyatnykh|19991|Zdravkov & ~ 


^ Zdravkov & Pamyatnylffi|200^ 

fference is the presence of lingers 


Saio|201 1 1. ITie notable difference is the presence of lingers 

in the lower boundaries of our j3 Cephei strips for ^ > 1. 
Their appearance here is due to the unprecedented resolution 
in HR-diagram space of our stability calculations. To eluci¬ 
date their origin. Figure plots part of the ( - \ frequency 
spectrum of an 8.5 Mq stellar model as it evolves from the 
ZAMS to the red edge of the main sequence (REMS), show¬ 
ing which modes are stable and which are unstable. The pi 
mode is unstable over the effective temperature range 4.358 > 
log(7’efr/K) > 4.317, and the gi mode over the cooler but over¬ 
lapping range 4.341 > log(reff/K) > 4.301. The star then 
passes through a phase with no unstable modes, before the in¬ 
stability reappears in the range 4.288 > log(7’eff/K) > 4.278 
for the g 2 mode. 

This alternation between instability and stability, seen as 
fingers in Figure!^ stems from the fact that the k mechanism 
only excites modes whose eigenfrequencies fall in a narrow 
range [ctiq, cry]. At frequencies Re(cr) > cr^i, the pulsation 
period becomes comparable to the local thermal timescale in 
the envelope region above the iron opacity peak, and this re¬ 
gion behaves as a damping zone, stabilizing the modes. Con¬ 
versely, at frequencies Re(cr) < cry, modes couple with grav¬ 
ity waves trapped in the p-gradient zone developing at the 
core boundary, and are likewise damped. The intermediate 
stable phase in Figure [T0| between log(7’eff/K) = 4.301 and 
log(7’eff/K) = 4.288 occurs when there are no modes in the 
[crjo, cry] range. As the star evolves, the unstable range nar¬ 
rows: cThi decreases due to lower T^ff, while cry increases due 
to the growth of the ju-gradient zone. 

Figure[^shows a version of the { - I panel calculated us¬ 
ing OP opacity tables rather than OPAL tables. There is an 
overall shift of the instability strips toward higher luminosi¬ 
ties, an effect already noted by jPamyatnykhj (j 1 999|l . The fin¬ 
gers persist with much the same structure, supporting the fact 
that they are physical effects rather than numerical artifacts. 
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Returning now to Figure the post-MS extension of the 
SPB strips has been attributed in the literature to features in 
the Brunt-Vaisala frequency which reflect gravity waves at 
the boundary of the helium core, preventing them from pen¬ 
etrating i nto the core and b eing dissipated by stron g radiative 
damping. |Saio et al.|(|2006|l and|Godit et al. (|2009|l argue that 
the necessary feature is an intermediate convection zone (ICZ) 
associated with the hydrogen-burning shell, but more recently 
[Daszyhska-Daszkiewicz et^ ( 2013| l have shown that even a 
local minimum in the Brunt-VSisala frequency is sufficient to 
reflect modes. In the present case, the empirical mass thresh¬ 
old M > 9 Mq required for formation of an ICZ coincides 
with the lower boundaries of the SPB strip extensions. In the 
lowest-mass models above this threshold, the ICZ vanishes 
shortly after its appearance, but it leaves behind a narrow re¬ 
gion with a steep molecular weight gradient. This gradient 
causes a spike in the Brunt-Vaisala frequency, which serves 
in a similar manner to prevent gravity waves from entering 
into the core and being dissipated. 

The correspondi ng post-MS extension of the fi Cephei s trips 
was first noted by |Dziembowski & Pamiatnykh| (|1993|l, but 
has not received much attention in the literature. Figure 
shows that this extension has a well defined lower boundary, 
much like the SPB stars although situated at slightly higher 
masses, M > 10.5 Mq. We have determined that the ex¬ 
tension is also a consequence of ICZ formation; the shift to 
higher masses arises because it appears that multiple convec¬ 
tion zones, rather than a single one, are necessary to reflect 
waves at the core boundary in the case of p Cephei pulsators. 

3.2. Asteroseismic Optimization 



Uosc mod 56.3 [/iHz] 


Figure 12. Echelle diagram for the subgiant star HD49385. Observed fre¬ 
quencies are shown as filled circles = 0), triangles (^=1) and squares 
{£ = 2); black horizontal lines indicate the Icr error bars. Calculated frequen¬ 
cies of the best-fit model are oveiplotted as the con'esponding open symbols. 


To illustrate the updated asteroseismic capabilities of MESA, 
Figure 


measured by |Deheuvels et al.| (|2010]l, and the correspond¬ 
ing frequencies of the best-fit model determined using the 
astero extension. The calculations follow the same proce¬ 
dure detailed in section 3.2 of Paper II; the only significant 
differences are that the initial mass, helium abundance, metal 
abundance and mixing length parameter are refined using the 
downhill simplex algorithm rather than the Hooke-Jeeves al¬ 
gorithm; oscillation frequencies are calculated using GYRE 
rather than ADIPLS; and the surface c orrections to frequen- 
cies are evaluated using equation 4 of [Ball & Gizon|P014| l 
rather than with the| |Kjeldsen et al. (2008 1 scheme. 

Comparing Figure |12| against ngure 8 of Paper II reveals 
only small differencesTetween the two. The of the best-fit 
models reported by astero is 2.3 in the former case, com¬ 
pared to 2.4 in the latter (cf. table 2 of Paper II). 

3.3. Automated Strip Location 




o 


12 plots the echelle diagram for the subgiant star 
HD 49385, showing both the frequencies of ^ = 0 - 2 modes 


l0g(reff/K) 


Figure 13. The growth timescale rose (left axis) and oscillation period Pose 
(right axis) of the fundamental and first-overtone radial modes of a 8.5 M© 
model, plotted as a function of Peff as the star evolves away from the ZAMS. 
The vertical dashed lines, determined automatically, show the points where 
the modes switch from stable (Tqsc < 0 ) to unstable (rose > 0)^ and vice versa. 

The instability strips presented above involved the exam¬ 
ination of ~ 11 million modes of ~ 40,000 stellar models. 
To partially automate the process we can leverage tight cou¬ 
pling between GYRE and MESAstar. This is achieved by mak¬ 
ing small modifications to the extras_checkjmodel callback 
routine in MESAstar (see Appendix B.l of Paper II), so that 
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Figure 14. The calculated blue edge of the classical instability strip, for fun¬ 
damental and first-overtone radia l modes. The correspondin g dashed lines 
show the predictions from set B of]Smolec & Moskalik|j2008| their figure 1). 


GYRE is run after each time step to determine the set of eigen- 
frequencies [cr] of a user-specified group of modes. When 
Im(cr) changes sign from one time step to the next for any of 
these modes, indicating that an instability-strip boundary has 
been crossed, a search is performed to find Im(cr) a; 0. 

Figure presents an application of the tight coupling 
to the fundamental and first-overton e ra dial modes of the 

8.5 Mq model considered in Section mi showing how the 
growth timescales Tqsc and oscillation periods Posc-l/vosc of 
the modes change as the star evolves from the ZAMS into 
the post-MS. The second-overtone radial mode remains sta¬ 
ble, Tosc < 0, over the range plotted. The vertical lines show 
where Tqsc changes sign. The blue and red edges of the {{ - 0) 
y6 Cephei instability strip can be seen in the upper panel of 
Figure [T^ at log(7’eff/K) = 4.36 and log(7’eff/K) = 4.30, re¬ 
spectively. The blue edgq^of the classical instability strip 
can likewise be seen in both panels at log(reff/K) = 3.75. 
The corresponding red edge is not found because GYRE does 
not include a treatment of the pulsation-convection interaction 
— a necessary ingredient for modeling the cl assical red edge 
(see, e.g, section 3.7.3 of |Aerts et al.||2OT0| and references 
therein). 

As a further demonstration of automated instability strip lo¬ 
cation, Figuref^plots the blue edges of the classical instabili¬ 
ty strip in the HR diagram, for fundamental and first-overtone 
radial modes. The edges are calculated for 51 evolution¬ 
ary tracks spanning the initial mass range 1.25 Mq < M < 

12.5 Mq in uniform logarithmic increments. At luminosi¬ 
ties log(L/Lo) > 2.5 corresponding to classical Cepheid pul- 
sators, these edges show good agreement with the set B results 
published by Smolec & Moskalik (|2008[ their figure 1). At lu- 


For classical (5) Cepheid pulsators, the observational blue edge of the 
classical instability strip is in fact established by stars evolving to higher 
Tefr on their first blue loop, rather than stars on their first crossing of the 
Hertzsprung gap. However, the purpose of the present section is to demon¬ 
strate the capability of tightly coupling, and in this context the distinction 
between the blue edges from multiple crossings is unimportant. 


minosities log(L/Lo) <1.6 corresponding to 6 Scuti stars, the 
edges are somewhat cooler than results published in the liter¬ 
ature; however, this is because we consider only fundamental 
and first-overtone modes, whereas the blue edge is typically 
set b y higher overtones w hich are displaced toward hotter Teff 
(e.g.,|Dupret et al.|2004| their figure 1). 

Ideally, the same automated approach could be used to lo¬ 
cate the boundaries of the non-radial {{ > 0) instability strips 
plotted in Figure[9] In practice it is very challenging to devise 
a robust algorithm that can unambiguously interpret the eigen- 
frequencies produced by GYRE. Sometimes, acoustic glitches 
in a model can trap modes in surface layers, where they are 
strongly excited; however, these modes are very sensitive to 
model parameters, and it is unclear whether they are physi¬ 
cally meaningful or not. 

4. IMPLICIT HYDRODYNAMICS 

Shocks happen in stars, such as after a massive star col¬ 
lapses, or cyclically in the outer envelopes of stars pulsating at 
sufficiently large amplitude. Previous versions of MESAstar 
allowed large velocities such as those encountered in the last 
few seconds leading to a core collapse (=» 1000 km s“'), but 
there was no provision for large jumps in velocities leading 
to shocks. In this section we describe the changes that have 
been made to support an implicit treatment of hydrodynamic 
shocks that includes careful attention to conservation of en¬ 
ergy. We demonstrate that the revised equations are intrin¬ 
sically conservative in the sense that deviations from exact 
energy balance can only arise from residual numerical er¬ 
rors in the approximate solutions rather than from the form 
of the equations themselves. Following the description of the 
changes, we show a series of envelope shocks as a test of the 
implementation. The form of the equations a nd the demon¬ 
stra tion of intrinsic co nservation closely follow Fraley ( |1968 1 
and Grott et al. P0051l. The tr eatment of artificial viscosity is 
based on |Weaver et al.| ( |1978| l. 
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Figure 15. Schematic of relevant cell and face variables relevant for hydro¬ 
dynamics in MESAstar. 


4.1. Mass Continuity 
The specific volume of cell k is 
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(25) 
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where is the outer face radius, r^+i is the inner face radius, 
dnik is the cell mass, and pk is the cell average density (see 
Figure 15 for the layout of cells in MESAstar). We create an 
initial algebraic form of the continuity differential equation 
by dividing the change in the specific volume in a step by the 
length of time 6t, using step start and end values for r^., 
and pk and the value dmi^ which is constant during the step; 


l/pk - l/Pstarg ^ 47r [(>-^ - rl^,k) " 

6t 3dmk 6t 

(26) 

Next, we rewrite the right hand side introducing variables 
for the time centered velocity v and the effective area A to 
get the final form of the mass continuity equation as used in 
MESAstar; 


more between neighboring cells. Making the shock spread 
directly dependent on the local cell widths would produce nu¬ 
merically intolerable dynamically changing, non-monotonic 
variations in cell-to-cell values for the shock spread. Use of a 
local ranning average cell width is also ruled out by the need 
to keep algebraic equations dependent only on nearest neigh¬ 
bors to allow a block tridiagonal matrix solution. The use of 
a small fraction of the local radius gives a smoothly varying 
shock spread that avoids the numerical problems associated 
with using the cell width. 

We define the quantity Qvisc./t (having dimensions of en¬ 
ergy), in cell k as 


id\isc,k ~ , I I 

dmk \rk r^-Fi/ 


(34) 


1/P<: l/Pstart.A: ^ r a A \ 

-- - 3 — Ak+iVk+\) , (27) 

of dnik 

where 

^k — i^k + t’start,t :)/2 (28) 

and rk is evaluated as 

rk = ^'start.f + Vkdt . (29) 

Algebraic simplification then shows that 

= y {rj + r^rstart,^ + rl^n,k) ■ (30) 

To be consistent with the mass continuity equation, we use 
these expressions for effective area and time centered veloc¬ 
ity in the following momentum and energy equations. It will 
be shown below that to get intrinsic energy conservation, we 
must time center the velocity and use special combinations 
of starting and ending radius in a couple of places, but all 
other terms in the equations can remain fully implicit to avoid 
degrading the numerical stability as would happen in a uni¬ 
formly time-centered scheme. 

4.2. Artificial Viscosity 

In MESAstar, the artificial dynamic viscosity coefficient 
77visc (which has the dimensions g cm^'s^') 

t]^[sc,k ~ ^visc,linear,t: t ^visc.quad.t: 5 (31) 

where the linear term is 


The momentum equation uses Qvisc in an expression that de¬ 
fines an artificial acceleration analogous to the pressure gradi¬ 
ent term, and the energy equation uses it to define an artificial 
viscous heating analogous to the mechanical work term. 


4.3. Specific Linear Momentum Equation 

The local linear momentum conservation equation at face k 
between inner cell k and outer cell k - 1 is 


^k r start,/: 

df 


GiUk 
fkf start,/: 


/-I 


- I ^visc,/ 

dmk 


(35) 


where dm^ - (dm^-\-dm}^^\)l2 is the mass associated with face 
k, and the viscous acceleration term at face k is 


^visc,^ 


/ 2visc,A:-l Q V 


rk 


dmk 


(36) 


The use of the product rj.rstart,/ in the denominator of the grav¬ 
itation term is necessary for intrinsic energy conservation as 
will be shown below. 


4.4. Specific Energy Equation 

The local energy conservation equation for cell k between 
outer face k and inner face k -H 1 is 

^/ ~ ^start,/ _ Lk — Lk+\ p (AkVk — Ak-AlVk+l\ 

6t dmk \ dmk j (37) 

^visc,/ ^nuc,/ ^v,k ^extra,/ 5 


^visc,linear,“ ^ ^lP//'mid,//'S,/ 


(32) 


where et is the specific thermal energy for cell k. The viscous 
heating rate for cell k is 


and the quadratic term is 

3 P^^did,/ 


'/vise,quad,^ ~ ^^2 


dnijc 


max 




r^+lV/+l 


■tivk 


) , (33) 


with rmidjt = {rk+i + rk)l2, c^^k the sound speed in cell k, and 
1 1 (I 2 ) is a dimensionless coefficient for the linear (quadratic) 
term. The linear term is rarely used; it provides for a general 
damping of pulsations. The quadratic term is only nonzero 
in regions of compression and is the primary control for the 
strength of artificial viscosity. Assuming the usual case of 
/i= 0, the shock fr ont is spread over a distance ~l 2 rk- We fol¬ 
low [Dorfi] (|1998)1 in opting for a shock spread proportional 
to the local radius r rather than the local cell width. This 
choice is dictated by the fact that step-by-step adjustments 
to the mesh resolution lead to dynamically changing, non¬ 
monotonic variations in cell widths of up to a factor of 2 or 


^visc,/ “ 


4^(2 vise,/ 

dmk 



V/+1 

^'z+i 


(38) 


Energy loss from weak reaction neutrinos is already sub¬ 
tracted from the nuclear burning term, fnuc,/, so only the neu¬ 
trino energy loss rate from thermal processes, Cy^k, is explicitly 
accounted for in Equation ([37|l. An example of fextra,/ is arti¬ 
ficial injection of energy to trigger a shock. 

An alternative form of the energy equation equates the 
model dLldm to the expected value 


where 


T/ Lk+i ^ IdL\ 

dmk /expected,/ 


(39) 


dL\ 

/ expected,^ 


“ ^grav,A: ^\isc,k ^nuc,A: ^v,k ^extra,/: 


(40) 
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and 


^grav,^: “ 

Using Equation ( |27l l, the expression for 
ek 


(41) 


^start.t: p l^k^k -^tr+lV^+l \ 

6t \ dnik ) 

can be rewritten 


-grav ’ 


^gra.y,k ~ 


^start,/: p (1 / Pk 1 / PstaiUk) 


6t 


5t 


(42) 


thereby avoiding the use of velocities and thus be appropriate 
for hydrostatic cases. 

4.5. Intrinsic Energy Conservation 
The summed kinetic, potential, internal energies are 


KE = ^ i dnik vl, 

k 

(43) 

pg, v Gnikdnik 

(44) 

IE = ^ e/: dnik , 

(45) 


and thus the total energy of the star is E = KE + PE + IE. 
We now explicitly demonstrate that the equations we solve 
are formulated in such a way that the rate of change of total 
energy exactly equals the combined energy sources and sinks. 

Multiplying Equation ( |T5| ) by Vkdifik gives an equation with 
units of luminosity: 


^ difik 


^ V? — 

k start, A: 

* 


= — Gnikdnik 

\fkf start,/; 

• AkVk{Pk-i - Pk) 


Vk 


(46) 




1 ^visc,/:)* 


Using Equation ( |29] l to eliminate Vk in the first term on the 
right. 


— Gnikdnik 


( V, ) 

1 Gnikdnik | 

fl 1 ] 

\^k^ start, ^ / 

^ St ' 

Vk ^ start, ^ / 


(47) 


shows that this term is the negative of the rate of change of po¬ 
tential energy, a result that is made possible by the use of the 
Gniklrkrstsii^k in Equation ( [rS]) instead of an alternative such 
as Gniklrf. Thus, Equation (|46[l can be written as 


1 -J- ''start,/: Gmkdnik I 1 1 

- dnik -I- 


6t 


- -AkVk{Pk-i 


6t 

Vk 


Vk 


^ Start,A: 


(48) 


' Pk) "t" 47r j ^ j (Qvisc,k-1 Qvisc,k) 


Similarly, multiplying equation ( [J7] i by dnik also yields an 
equation with units of luminosity: 


i^k ^start,/:) 

Jt 


dnik =-iLk- Lk+i) 

- Pk (AkVk - Ak+\Vk+\) 

. „ (Vk V/t+l 

+ 47r^visc,A: I 

Vk 

(^nuc,^ “ G,k ^extra,fe) dnik 


(49) 


Adding Equations ( |48] l and ( |49] l and summing over the grid 
index k gives 


E Gnikdnik (I 1 

6t \fk Atart,/; 


1 T- 'i - ‘'L.i 

*2 -^r- 

{Ck ~ ^start,A:) , 

+-;- dnik 

ot 


- ~ Pk+\) 


(50) 


+PkAk+\Vk+\ -Pk-iAkVk 

+47r2visc,,t-i I — 1 - 4;r2visc,4 — 

\rk) \rk+i 

"l"(^nuc,A: “ ^v,k ^extra,/:) dlflk ■ 

The sum over the pressure terms is 

(Pk^+iVk+i - Pk-l^kVk) 
k 

~ ~ [(7^^i^)surface “ (7^^i^)center] 

— (Tacoustic,surface Tacoustic,center) ? 


(51) 


which cancels term by term except at the boundaries. We de¬ 
fine Tacoustic.surtace as the work done by the model on the at¬ 
mosphere at the surface and Tacoustic.center as the work done on 
the model at the center, for example, by an artificial piston. 
The sum over the artificial viscosity terms leads to a similar 
expression, but because Qvisc vanishes at the surface and the 
center, the sum equals zero. That is, the energy added by arti¬ 
ficial viscous heating in the energy equation exactly balances 
the loss of kinetic energy by artificial viscous acceleration in 
the momentum equation. 

The terms on the left hand side of Equation ( [50l l are the dif¬ 
ference in the total energy between the start and end of a step 
divided by the length of the step, in other words, the average 
rate of change of the total energy of the model. Therefore 
Equation ( |50| can be written as 

(Efinal “ E'lYiiii^l) /St — — (Usurface Teenier) 

“ (Tacoustic,surface Tacoustic.center) 

+ ^ (fnuc,/: - ^v,k + ^extra,/:) dnik . 
k 

(52) 

This equation embodies conservation of energy in MESAstar: 
the rate of change of total energy equals the combined en¬ 
ergy sources and sinks. This demonstrates that in the given 
form, the algebraic equations intrinsically conserve energy 
in the sense that failure to get energy balance can only arise 
from the residual numerical errors that are inherent in using 
approximate solutions to the equations. This in turn means 
that to control energy balance errors, we can focus on reduc¬ 
ing residuals either by changes in the Newton solver or by 
timestep reductions. 

4.6. Controlling the Accuracy of Energy Conservation 
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The Newton solver considers both the sizes of incremental 
changes to the variables and the sizes of residual errors for 
the equations. For the energy equation, the residual used by 
the solver is defined to be the timestep 6t times the difference 
between the left and right sides of Equation ( [JT) ! divided by 
^start,*:; in Other words, the residual is the error as a fraction 
of the specific energy at the start of the step. By adjusting 
tolerances for the average and maximum size of residuals, we 
force the Newton solver to take extra iterations to reduce the 
residuals which will in turn reduce the total error in energy 
conservation. 

A second and related way to control energy conservation 
errors is to use the average and maximum energy residuals 
to adjust the timesteps. For example, if the maximum mag¬ 
nitude for an energy residual exceeds a specified hard limit, 
then the proposed solution is rejected and the step is retried 
with a smaller timestep. If the maximum is smaller than the 
hard limit but exceeds another specified limit, then there is no 
forced retry for this step, but the next timestep is reduced by 
the ratio of the limit divided by the maximum magnitude. If 
the maximum is smaller than both limits, then other factors 
determine the next timestep. Later in this section, we will 
show that these approaches in combination with the intrinsic 
conservation form of the equations yield a solution for a shock 
in an envelope that evolves with reasonably large timesteps 
while conserving energy to a high degree of accuracy. 


4.7. Limiting Acceleration of Convective Velocity 

When using hydrodynamics, we often require timesteps 
that are so small that we need to limit the increase in convec¬ 
tive velocities as calculated in the standard instantaneous mix¬ 
ing length theory (MET) so that they do not assume unphysi- 
cally large accelerations. If convection velocities are allowed 
to adjust instantaneously, then our methods for artificially cre¬ 
ating shocks will fail since however rapidly we inject energy 
over a limited region, convection will be able to transport the 
energy away. To be able to simulate shocks we need to have a 

way to limit convection velocity acceleration. _ 

The prim ary scheme we u se for this is derived from |Ar-| 
|nett| ( |1969| l and [Wood] \\91A) . The MET implementation in 
MESA has been extended to take as additional arguments the 
timestep and the previous convection velocity at the same 
mass location (Vc.prev)- It calculates a provisional convection 
velocity (vco) using the standard instantaneous MET, then de¬ 
fines a convective timescale (tmlt) as the local pressure scale 
height (H) divided by the sum of the provisional plus previ¬ 
ous velocities. If 6t is less than tmlt, then the next convection 
velocity (v^) is only incremented from the previous one by 
the difference of the provisional minus the previous velocities 
times the ratio of the timestep divided by the time scale 

■ I \ 

— ^c.prev + min I 1, I (VcO “ Vc^prev) 5 (53) 

\ Tmlt/ 


where 


Tmlt - 


H 

(VcO + Vc.prev) 


(54) 


As an alternative scheme for limiting convection acceleration, 
we also allow the maximum rate of change of convection ve¬ 
locity to be set as a fraction, gg, of the local gravitational ac¬ 
celeration. If Vco > Vc.prev, then 

(55) 


The final Vc is used to recalculate the convection efficiency, 
which is used to calculate the MET temperature gradient. 

These methods for limiting the acceleration of convective 
velocities reduce the energy transport rate as well as the rate 
of compositional mixing. Both schemes seem to give at least 
qualitatively reasonable results and avoid the problems of un- 
physically large accelerations that are possible with standard 
instantaneous MET. Hopefully this ad hoc solution will soon 
be replaced by a quantitatively correct formulation. 

4.8. Surface Boundary Conditions 

MESA provides a variety of options for surface boundary 
conditions (see, e.g., section 5.3 of Paper I), and several more 
have been added for use with hydrodynamics. The simplest 
allow specification of a particular value for the surface pres¬ 
sure, the surface temperature, or the Tetr if the surface is not 
at the photosphere. In the case of a given fixed surface pres¬ 
sure, the corresponding surface temperature is set using the 
surface luminosity and radius based on the usual black body 
relation. For the second case, where the surface temperature 
is fixed, the surface pressure is set to the corresponding ra¬ 
diation pressure. For both of these, if the surface is not at 
the photosphere, Teff is set using the Eddington T -t relation. 
Finally, for specified Teff when the surface is not at the photo¬ 
sphere, the corresponding surface temperature is also derived 
using the Eddington T -r relation, and the surface pressure is 
set to the radiation pressure for that temperature. 

For computations involving shocks at the surface, there is 
an option to use boundary conditions that specify a vanish¬ 
ing gradient for compression at the surface and a temperature 
corresponding to black body radiation. The outermost cells 
{k - 1,2) satisfy the equation 


Pi Pstart.l P2 Pstart,2 

which represents the vanishing of the surface compression 
gradient. 

Finally, for computations involving interior shocks but low 
velocities at the surface, there is an option to use the surface 
pressure from the selected atmosphere prescription with the 
momentum equation relating the surface velocity to the sur¬ 
face pressure gradient. This form for the surface boundary 
conditions is used in the shocked massive star example in Sec- 
tionj^and in the following envelope shock test. 

4.9. Shock Test 

To test the implementation, we shock the extended enve¬ 
lope of a 6.93 Mq asymptotic giant branch (AGB) star evolved 
from a 7 Mq main-sequence star without rotation and an ini¬ 
tial metallicity of 0.0()1. This case is chosen because of the 
uniform properties of the extended envelope (i.e., small den¬ 
sity range, smooth density, and uniform composition). Our 
interest is to study the propagation of the shock, the proper¬ 
ties of the shocked material, and the magnitude of energy con¬ 
servation errors. In Section]^ we present results that mimic 
core-collapse supemovae. 

Explosion simulations with MESA start from a converged 
model. The core is excised by removing inner shells of the 
model and setting new inner boundary conditions for mass, ra¬ 
dius, velocity, and luminosity. For the current test, we remove 
the center just above the helium core at a mass of 2.40 Mq 
which corresponds to an inner radius of 27.2 Rq. The stel¬ 
lar surface lies at a radius of 282.7 Rq. During the following 


Vc = min(VcO, Vc.prev + dt gg g) , 
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Figure 16. Multi-epoch snapshots for the hydrodynamical simulation of a 10^^ erg shock in the envelope of a 6.93Mo AGB star. We show the density (top row), 
temperature (middle row), and velocity (bottom row), versus Lagrangian mass (left column) and radius (right column). In each panel, the solid line refers to the 
MESA results and the dashed line to the VID results. 


evolution, the excised region is treated as a point mass and is 
linked to the above layers by the inner boundary conditions 
which can be changed at each timestep to simulate various 
core events. The model grid was adjusted at each step to give 
higher resolution in the vicinity of the shock. The total num¬ 
ber of cells stayed at about 1000, with cell masses dropping 
to about 10“^ of the total in approximately 100 cells around 
the shock. 

In MESA, the artificial explosion that creates a shock can 
be produced in three ways: a piston, a luminosity flash, or a 


thermal bomb: 


The first option changes the inner boundary conditions 
for velocity and radius to mimic a piston. A core¬ 
collapse supernova can be simulated by moving the in¬ 
ner radius inwards (collapse) at a free-fall speed and 
then violently outwards (bounce and explosion). The 
parameterization for the p iston-driven explosion is th e 
same as that described in|Woosley & Weaver] (|1995|), 
and includes the infall piston time, the hnal inwi'd pis- 
ton radius, the initial outward piston speed, and the final 
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piston radius. 

• The second option increases the inner boundary lumi¬ 
nosity over a specified time in order to deliver the de¬ 
sired total energy. In this approach, the inner boundary 
radius is fixed at all times and becomes a zero-flux in¬ 
ner boundary once the flash is over. In this approach, 
the inner boundary radius is fixed (zero velocity) and 
we inject the energy within the first zone of the domain. 

• The third option deposits energy at a constant rate dur¬ 
ing a specified time and in a region bounded by two 
specified Lagrangian-mass coordinates. As in the sec¬ 
ond option, the inner boundary radius is fixed at all 
times. 

The differences amongst these three options can alter the 
properties of the shocked envelope. 

To benchmark MESA for these shock tests, we have used 
the explicit radiation-h ydrodynamics code VID ( |Livne|[T993[ 
[Dessart et al.|2010a|b| l. Options 1 and 3 are implemented in 
VID. For the present envelope shock test, and subsequently 
for the explosion tests, we use option 3 in both codes. We 
initiate the explosion by depositing a total of 10"^® erg at a 
constant rate over 10 s between the Lagrangian mass coordi¬ 
nates of 2.40 and 2.45 Mq. This energy deposited is well in 
excess of the initial binding energy, which is approximately 
-2x10^^ erg. Once the energy injection is over, we save a 
model which is then used as initial conditions for a shock evo¬ 
lution simulation. 

Once the stellar core has been excised, the remaining en¬ 
velope has a smooth density profile, resembling a power law 
whose exponent is -1 at depth and decreases outwards to be¬ 
come about -10 at the photosphere (top row panels of Fig¬ 
ure [^. Because convective accelerations are limited, the en¬ 
ergy deposited increases the internal energy within the inner¬ 
most 0.05 Mq of the grid. The pressure build-up leads to the 
sudden expansion of the innermost layers and the formation 
of a mildly supersonic shock (Mach number a; 2). The shock 
propagates at a velocity in excess of 1000 km s“' initially, but 
slows to a few 100 km s“* by the time it reaches the stellar sur¬ 
face after 3x10^ s. The density contrast across this somewhat 
weak shock is ^ 6. For a strong shock, one expects a density 
jump of 4 for an ideal gas with an adiabatic index of 5/3 and 
a value of 7 for a radiation-dominated gas (y - 4/3). 

This simulation is analogous to a shock-tube test. However, 
in the stellar context (realistic stellar envelope, realistic equa¬ 
tion of state, spherical expansion), there is no analytical so¬ 
lution for comparison. We thus run the same simulation with 
the code VID and include the results in Figure [T6| The results 
agree at multiple times spanning the progression of the shock 
towards the stellar surface (the times used for comparison are 
the same to within 1% and the grid resolution is comparable). 
The sharpness of the shock in the two simulations differs with 
time and location. In VID, the artificial viscosity has a phys¬ 
ical spread of two grid zones, irrespective of radius, while in 
this ME SA ru n, the spread is set to 0.1% of the local radius (see 
Section[4^. 

Since the explosion is started as a thermal bomb, the bulk 
of the energy is initially internal, see Figure [T^ As the ma¬ 
terial expands and accelerates, the kinetic energy increases, 
mirroring the decrease in internal energy (essentially no en¬ 
ergy is used to unbind the envelope). At the time of shock 
emergence, the internal and kinetic energies are comparable. 



t [s] 


Figure 17. Top: Evolution of internal energy Eint, gravitational energy Egrav, 
kinetic energy Ekin. and their sum Etot for the envelope shock test simulated 
with MESA and VID. Bottom: Log of cumulative relative error CRE (Equa¬ 
tion of the total energy Etot- We neglect sources (nuclear burning) and 
sinks (radiation losses), which are negligible. 
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Figure 18. Top: Normalized values of velocity, artificial acceleration gvisc. 
artificial viscous heating 6visc. and Mach number in the vicinity of the shock. 
The dashed vertical line marks where the Mach number is unity. Bottom: 
Dependency of the shock morphology on changes in the viscosity parameter 
I 2 . The dots shown for the model with I 2 =0.004 denote the location of the 
MESA grid points at that time. For all these tests, we deposit an energy of 
10 ^^ erg at a constant power over 1 s. 


In the present case, we can preserve good accuracy while 
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still allowing time steps an order of magnitude greater than the 
Courant time{^ The error in energy conservation at timestep 
i is 

dfi” error,i “ ^/-l ^sources,/ (57) 


where fisources,; is the right hand side of Equation ( |52] i mul¬ 
tiplied by St. The cumulative relative error in energy at a 
timestep n is 


1 " 

CRE{t„) = — YsE, 


F ^ 

l=l 


(58) 


In the test case, after about 15,100 timesteps when the shock 
reaches 6.6 Mq, the cumulative relative error has grown to 
about -1.4 X 10^®, corresponding a roughly linear growth rate 
of about -1X10“'° per step (bottom panel in Figure 17 1 . Note 
that at this stage of evolution the shock is nearing tEe outer 
edge of the envelope but has not actually broken out through 
the surface. Issues of shock break out are beyond the scope 
of the current implementation. Using the parameters selected 
for the test, the energy conservation with VID is not as good as 
with MESA (the jumps in cumulative error correspond to times 
when the limit on the time step are loosened); comparable 
accuracy can be obtained with VID by reducing the explicit 
time step well below the Courant limit. 

Finally, to illustrate the effects of artificial viscosity, we 
vary the quadratic term I 2 that controls the spread of the shock 
in response to compression (see Equation [3^, with the ex¬ 
plosion energy increased to 10^° erg in order to produce a 
stronger shock, and otherwise the same parameters and ini¬ 
tial conditions. The top panel of Figure 18 shows the artifi¬ 
cial acceleration (gvisc) and energy (fvisc) terms that enter the 
momentum and the energy equations for I 2 =0.001. The ac¬ 
celeration term is positive ahead of the shock, causing a pre¬ 
acceleration of the unshocked material, and negative behind 
the shock causing a deceleration of the post-shock material. 
The energy corresponding to those changes in momentum is 
balanced by the extra term for artificial viscous heating in the 
energy equation (evisc). The lower panel of Figure 18 shows 


the expected increase in the width of the shock as we raise 
the parameter l 2 . For the model with I 2 =0.004, dots locate 
grid cells. Note that with the smallest value {I 2 =0.001), the 
velocity is showing small oscillations (“ringing”) behind the 
shock indicating that we have reached a practical lower bound 
for the shock spread given the other parameter choices and the 
nature of the specific problem. 

5. ADVANCED BURNING 

For the advanced stages of stellar burning, we show here 
that more accurate summations yield more efficient time in¬ 
tegrations. This development allows MESA to use large in- 
situ reaction networks. It offers an improvement by providing 
a single solution methodology that avoids the challenges of 
stitching together different solution methods such as nuclear 
statistical equilibrium (NSE) or co-processing a reaction net¬ 
work. We discuss this development and apply it to the evo¬ 
lution of an X-ray burst on a neutron star. In Section we 
discuss the pre-supernova progenitors and combine the new 
capability for advanced burning with the implicit hydrody¬ 
namics module to discuss the explosion of core-collapse su¬ 
pernovae. 

The Courant time, equal to the minimum sound crossing time through a 
grid zone. In this envelope test, it is of the order of 10 s initially, increasing 
progressively to 40 s prior to shock emergence. 


The equations that describe the continuum limit of reacting 
nuclei are 




J 

+ 


(59) 


where T, is the abundance of isotope i, 4 is a reaction rate, and 
the three sums are over reactions which produce or destroy a 


nucleus of s 


recies i with 1, 2, and 3 reacting nuclei, respec- 


tively (e.g.. 

Meyer et al. 19981 Hix & Meyer 2006( Guidryj 

et al.||20131 

Longland et al.||2014i. ITie positive or negative 


stoicniometnc coemcients c, account for the numbers of nu¬ 
clei created or destroyed in a reaction. The factorials in the 
denominators avoid double counting of identical particles. 


1 ITT I I I I I I I I I I I I I I I I T 
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Figure 19. Evolution of the composition for a one-zone burn at constant 
T=9.6xl0^ K and p=6.0xl0^ g cm“^ for 10^ s starting with a pure ^^Si 
composition. The calculation uses the mesa_204.net isotope listing (see 
Tablejg, the most abundant isotopes are drawn with thick lines, and several 
isotopes are labeled. The initial composition is quickly erased as NSE for 
Ye » 0.5 is established by « 10“^ s. Several orders of magnitude in time pass 
before weak reactions drive a second period of reaiTangement. By % 10 s a 
second NSE quiescent period with Yq % 0.403 is established. 


Figure [T^ shows the evolution of the mass fractions for a 
MESA one-zone burn at constant T =9.6x10° K and p=6.0xl0° 
g cm“° for 10° s starting with a pure ^^Si composition. The 
204 isotope network, mesa_204.net, used in the calculation 
is listed in Table 7\ and includes the isotopes identified in 
|Heger et al.| ( |2001 i as important for in core-collapse mod- 
els. The thermonuclear re action rates are from JINA reaclib 
version V2.0 2013-04-02 ( |Cyburt et al.||2010) l. Implementa¬ 
tion of reaction rates and associated quantities are described 
in Paper I and Paper II. 

The thermodynamic conditions used in Figure 
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are 

representative of the central regions of massive stars dur¬ 
ing the advanced stages of evolution. At such temper¬ 
atures the initial composition of pure ^^Si undergoes a 
rapid readjustment. The timescale for the initial Yg « 0.5 
composition to relax to an NSE composition is roughly 


exp(179.7/r9 - 40.5) s = 3 x 10“'' s (Khokhlov 


TO^ICalderetaT 


dJ2007| 
191 Bet' 


phase in Figure 
topes "'He and dominate the U 


, commensurate with the first burning 
etween « 10“^ s and =» 10“"' s the iso- 
0.5 NSE composition. 
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Table 2 

204 Isotope Network Listing. 


Element 

■^min 

^max 

Element 

^miii 

^max 

n 



s 

31 

37 

H 

1 

2 

Cl 

35 

38 

He 

3 

4 

Ar 

35 

41 

Li 

6 

7 

K 

39 

44 

Be 

7 

10 

Ca 

39 

49 

B 

8 

11 

Sc 

43 

51 

C 

12 

13 

Ti 

43 

54 

N 

13 

16 

V 

47 

56 

O 

15 

19 

Cr 

47 

58 

F 

17 

20 

Mn 

51 

59 

Ne 

19 

23 

Fe 

51 

66 

Na 

21 

24 

Co 

55 

67 

Mg 

23 

27 

Ni 

55 

68 

A1 

25 

28 

Cu 

59 

66 

Si 

27 

33 

Zn 

59 

66 

P 

30 

34 





Table 3 

Final Fe for Figure [T^ 


# of Isotopes 

U 

^max 

^max 

75 

0.4093 

Ni 

68 

125 

0.4065 

Ni 

68 

160 

0.4032 

Ni 

68 

204 

0.4032 

Zn 

66 

368 

0.4035 

Zn 

77 

833 

0.4029 

Sn 

125 

3298 

0.4039 

At 

211 


Since T and p are constant, only changes to Y^. can change the 
abundances. A second period of intense rearrangement begins 
at ai 10“"^ s and ends at a; 10 s. This activity is driven primar¬ 
ily by p{e^,v)n and n{e^,v)p and other weak reactions that 
change Y^. Beyond a; 10 s the isotopes '^^Ca, '^^Ca, and ^'Sc 
dominate the Fe ~ 0.403 NSE composition. 

Table shows the sensitivity of the final Fg in this calcula¬ 
tion to the number of isotopes in the network. Each succes¬ 
sively larger network encompasses the previous smaller net¬ 
work and was crafted to yield approximately the same final 
Fe value as given by the largest network. The 204 isotope net¬ 
work used in Figure[^is in the regime where larger networks 
give the same final Fg to 3 significant figures. 

5.1. More Accurate Summations Yield More Efficient 
Integrations 

We now test different summation methods for Equation ( [59] l 
and demonstrate that improved accuracy of the summations 
reduces the number of time steps with a commensurate reduc¬ 
tion in the execution time - while producing the same answers 
to within the specified integration accuracy. 

When the summations in Equation ( [59| are accumulated in 
IEEE 64-bit arithmetic (16 significanthgures, real*8 pre¬ 
cision in Eortran on most architectures; more specifically bi- 
nary64 with round to nearest and round ties to even) the in¬ 
tegration in Eigure [Intakes 3062 time steps using a variable- 
order Bader-DeuflhSd integrator with a specified accuracy of 
Tint=10“^ and a scaling value yscaie=10“ • The specified ac¬ 
curacy Tint limits the maximum error over one time step for 
any isotope. Other potential, but less demanding, choices 
for the meaning of Tint include limiting the average or root- 
mean-square error over one time step for all isotopes. When 


an abundance is greater than yscaie a relative error is calcu¬ 
lated, while for abunda nces smaller than yscaie. the absolute 
error is calculated (e.g.. Press et ar||1992| l. In essence, only 
abundances greater than yscaie can exert control on the size of 
the time step. 

When the summations are accumulated in IEEE 128-bit 
arithmetic (32 significant figures, real* 16 precision in Eor¬ 
tran on most architectures), the same integration takes only 
55 time steps, a factor ^50 improvement in the number of 
time steps, and a factor of =s30 less execution time. Both cal¬ 
culations returned the same answers to within the specified 
integration error tolerances. Eor tighter integration tolerances 
of Tint=10“® and yscaie=10“®, the evolution with summations 
in IEEE 64-bit arithmetic takes 10,081 time steps while the 
evolution with summations in IEEE 128-bit arithmetic takes 
88 time steps. This is a factor of a; 100 improvement in the 
number of time steps, a factor of a; 150 in execution time, 
with both calculations again producing the same abundances 
to within the specified integration error tolerances. Both sets 
of integration tolerances are practical, everyday usage toler¬ 
ances; they are not extreme cases of hypothetical interest only. 
Using low-order Rosenbrock and first-order Euler integrators 
also showed similar improvements in the number of time steps 
when the summations were performed in IEEE 128-bit arith¬ 
metic instead of IEEE 64-bit arithmetic. We achieve a reduc¬ 
tion in the number of timesteps and execution times regardless 
of the number of isotopes, choice of integrator, integration 
tolerances, or linear algebra solver. This improvement in effi¬ 
ciency is fundamentally driven by a reduction in the numerical 
noise of the function being integrated. 

At temperatures larger than 5x10® K, integrating Equa¬ 
tion (591 can be challenging as terms in the summation usu¬ 


ally become large and opposite in sign. As shown above, the 
classic symptom during an integration under these thermody¬ 
namic conditions is the integrator taking an excessive number 
of very small time steps in order to satisfy the specified inte¬ 
gration accuracy criteria. The traditional workaround to this 
numerical problem is abandoning a network integration at ele¬ 
vated temperatures and deploying equilibrium solution meth¬ 
ods. This switching of methods raises its own numerical is¬ 
sues when used within the larger context of multi-di mens ional 
simulations or stellar evolution models (see Section |5^ . 

Unless precautions are taken the summation of large sets 
of numbers can be very inaccurate due to the accumulation 
of rounding errors. Methods for accurate summation within 


search (e.g., Demmel & Hida 

20031 

McNamee 20041 Ogita 

et al. 

2005) Rump et al.|2008 

Grail] 

at & Menissier-Morain 

7un 

Collange et al.|20I4|l. ITiese summation discrepancies 


also worsen on heterogeneous architectures - such as clus¬ 
ters with NVIDIA GPUs or Xeon Phi accelerators - which 
combine programming environments that may obey various 
floating-point models and offer different precision results. 

The summations in Equation ( [59] l for the neutron, proton, 
and a-particle abundances are especially prone to inaccura¬ 
cies because every isotope in a network reacts with these three 
particles. We report on the summation errors for these three 
isotopes. Each term in the summations of Equation (|59]l is 
calculated using IEEE 64-bit arithmetic and then copieointo 
a IEEE 128-bit variable using the Eortran promotion rules. 
Each IEEE 128-bit term is then imported into the MP (|Brent| 
1978|l and MP£90 (|Bailey| 1995)1 multiple precision packages. 
All other aspects of the integration were executed in IEEE 64- 
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Table 4 

Results of Summation Experiments. 


IEEE 

Arithmetic 

Maximum 
Digits Compared 

Strategy 

Minimum 
Con'ect Digits ^ 

Number of 
Timesteps^ 

Ratio of 
CPU Times" 

ri„t = 10 Vscale = 10 ^ 

64-bit 

16 

in order given 

6 

3062 

31.7 

64-bit 

16 

sorted, ascending 

7 

2614 

24.4 

64-bit 

16 

sorted, Kahan sum 

8 

1141 

13.1 

128-bit 

32 

in order given 

21 

55 

1.0 

128-bit 

32 

sorted, ascending 

22 

55 

1.0 

rint = 10 ^ .Vscale = 10 ^ 

64-bit 

16 

in order given 

6 

10081 

156 

64-bit 

16 

sorted, ascending 

7 

7972 

123 

64-bit 

16 

sorted, Kahan sum 

8 

7674 

112 

128-bit 

32 

in order given 

21 

88 

1.0 

128-bit 

32 

sorted, ascending 

22 

88 

1.0 


^ Relative to the 100 digit sum by the MP and MP£98 multiple precision packages. 

** For a Bader-Deuflard integrator in IEEE 64-bit arithmetic. 

For a single thread on one 2.7 GHz Intel Xeon E5 core with the Intel 15.0.1 Fortran compiler, and relative to the 
execution time for the integration with 128-bit summations with terms in the order given. 


bit arithmetic. The summations are then accumulated with 

• IEEE 64-bit; terms in the order as given. 

• IEEE 64-bit; terms sorted by their absolute value in as¬ 
cending order. 

• IEEE 64-bit; terms so rted by their ab solute value in as¬ 
cending order and the |Kahan| (|1965|l algorithm, which 
reduces the numerical error in summation by retaining 
a separate variable to accumulate the errors. 

• IEEE 128-bit; terms in the order as given. 

• IEEE 128-bit; terms sorted by their absolute value in 
ascending order. 

• MP and MP£90 100 digits; terms sorted by their absolute 
value in ascending order. 

There are many summation methods and alternative multiple 
precision packages we did not deploy in these studies (e.g., 


Knuth|19971|Higham|2002HLi et alpOOH [Muller et al.|2010t 

Collange et al.|2014|l. ' 

Table|^summarizes these summation experiments. We con¬ 
firm that 100 digits are sufficient to prevent errors in our mul¬ 
tiple precision sums. Column 4 gives the minimum number 
of correct digits in a summation. There are three time periods 
in the evolution of Eigure where summations performed 
in IEEE 64-bit arithmetic greatly increase the number of time 
steps taken by the integration. One is during the first rear¬ 
rangement into the NSE state ending around 10“^ s, another is 
during the second rearrangement around 10 ' s and the third 
time period is when the abundances do not change much (T =» 
0) and reaction rates reach equilibrium. It is during these equi¬ 
librium periods where a summation in IEEE 64-bit arithmetic 
with the terms summed in the order given may yield only 6 ac¬ 
curate digits. As a result, 3062 time steps are needed to com¬ 
plete the integration (e.g., row 2 of Table It is important to 
note that this strategy and choice of ari thmetic is com monly 

rimme s|1999 1 


used by nuclear reaction networks (e.g.. 


and 


the most inaccurate choice. Row 4 of Table is an important 
case, sorted plus Kahan summation, because it demonstrates 
that a marginal improvement in the accuracy of the summa¬ 
tion (8 minimum correct digits) has a major reduction on the 
number of time steps (1174 time steps) and execution time (a 
factor of ss2.5 smaller). This establishes the general trend that 
improved accuracy of the summations reduces the number of 
time steps with a commensurate reduction in the execution 


time - while producing the same answers to within the spec¬ 
ified integration accuracy. 

The left hand side of Equation (59 1 for T,- is a IEEE 64-bit 
array to be filled with one of the summations. Setting T, equal 
to one of the IEEE 128-bit summations (at least 22 digits of 
accuracy) or the multiple precision package summations gives 
the most efficient integration (55 time steps, rows 5 and 6 in 
Table because the Eoitran precision demotion rules assure 
Yi is accurate to the limit of IEEE 64-bit arithmetic. The next 
best strategy, but a distant second, is setting T, equal to the 
sorted, Kahan summation. The worst case is setting T, equal 
to the 64-bit arithmetic sum with the terms in the order the y 
appear - which is a common approach (e.g., Timmes|1999|l. 



Figure 20. Number of accurate digits in 64-bit and 128-bit summations for 
y('^He) as measured by the 100 digit sum calculated by the multiple precision 
packages. The x-axis gives the time step number for the integration done with 
IEEE 128-bit summations. The number of accurate digits in y(p) and y(n) 
are within a few digits of yf^He). 


Figure shows the number of correct digits in 64-bit and 
128-bit summations for T("^He) with the terms accumulated 
in the order they are given. The number of correct digits is 
measured against the 100 digit sum calculated by the multiple 
precision packages MP and MPf9Q. The choices for the inte- 
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grator and integration tolerances are the same as in Figure [T9| 
Figure [ 20 ] shows that the minimum number of accurate digits 
is usually within a few digits of the limit of IEEE 64-bit arith¬ 
metic, but degrades to 6 digits (see row 1 of Table ^ during 
a time period of intense isotope rearrangement. These rela¬ 
tively large inaccurate summations cause the right hand side 
of Equation ( [59| to be poorly defined in IEEE 64-bit arith¬ 
metic. As a direct result, the integration of Equation ( [59| with 
IEEE 64-bit summations takes 3062 time steps to complete. 
Sorting the terms in the sums in ascending order and using 
the Kahan summation algorithm results in an accurate digit 
pattern that is very similar except the number of accurate dig¬ 
its is improved by one or two (see row 3 of Table |^. As a 
direct result of the improved accuracy of the summations, the 
number of timesteps is reduced from 3062 to 1141. 

Eor the IEEE 128-bit sum relative to the 100 digit sum in 
Eigure the minimum number of accurate digits is usually 
near thehmit of IEEE 128-bit arithmetic, but degrades to 21 
digits (see row 5 of Table during the second period of in¬ 
tense isotope rearrangement. Relative to the IEEE 64-bit sum¬ 
mations the number of correct digits is improved by at least 
15, consistent with our conversion of the IEEE 64-bit terms 
in the sum to IEEE 128-bit using the Eortran promotion rules. 
As a direct result of the improved accuracy of the summa¬ 
tions the integration that used the IEEE 128-bit summation 
took only 55 timesteps to complete the evolution. 

We found no notable improvements by increasing the accu¬ 
racy of the summations for the Jacobian matrix used by the 
stiff ordinary differential equation integrators. Eor this prob¬ 
lem, it is evidently more important to better define the func¬ 
tion - the right hand side of Equation ( [59] l - than the Jacobian 
matrix holding the derivatives of the function. 

Based on these experiments we currently choose to improve 
the accuracy of the summations by converting the terms in 
the order they appear from IEEE 64-bit to IEEE 128-bit us¬ 
ing the Eortran promotion mles and adding the terms in IEEE 
128-bit arithmetic. IEEE 128-bit precision is presently almost 
always implemented in software by a variety of techniques 
(e.g., double-double methods), since direct hardware support 
for IEEE 128-bit precision is presently rare. However, Table 

shows the reduction in the number of time steps from accu¬ 
mulating the sums in IEEE 128-bit arithmetic far exceeds the 
extra computational cost per addition. 

5.2. A Uniform Solution Methx)d for Nuclear Burning in 
Stellar Evolution 

At high temperatures, the traditional workaround for the 
numerical problem of inaccurate summations in IEEE 64-bit 
arithmetic is to forgo using a reaction network integration to 


above 3x10® K, and then move into NSE range above 5x10® 
K. Ad-hoc decision trees must be created for switching be¬ 
tween a network integration, QSE solutions, and NSE solu¬ 
tions. These switches can introduce unphysical discontinu¬ 
ities in the abundances either from one timestep to the next or 
in the abundance spatial profiles from one cell to the next. 

Eurthermore, cells near the transition between a network 
integration and an equilibrium method can be unstable in the 
sense that the equilibrium solution can evolve a cell to lower 
temperatures pushing the cell into using a network integration, 
while the solution from the network integration can evolve 
the cell towards higher temperatures evolving the cell back 
towards using the equilibrium solution. Moreover, the reac¬ 
tion network used for the time integration is different (usu¬ 
ally smaller) than the isotope listing used for the equilibrium 
methods. This necessitates crafting a delicate mapping be¬ 
tween two abundance vectors, which may also introduce un¬ 
physical discontinuities. In addition, care must be taken to 
assure the reaction rate screening corrections used in the time 
integration are properly taken into account in the equilibrium 
solution method, otherwise a fundamental incompatibility ex¬ 
ists between the abundance vectors. 

Einally, equilibrium methods determine the composition at 
a fixed electron fraction Tg- It then becomes necessary to 
solve an ordinary differential equation for Te based on weak 
reaction rates in or der to advance the abundance solution with 
a time varying Tg (McLaughlin et al. 1996t [Townsley et al. 


|2009[|Arcones et al.|2010| also see Section|8|. Switching be 

tween integration and equilibrium methods mid-stream is a 
liability, not a positive asset. 

The need for traditional workarounds forced by limited ac¬ 
curacy of the summ ation is now avoided. The summation ex¬ 
periments in Section 5.1 demonstrate that network integration 
can be robust and efficient, even at very high temperatures, 
when the accuracy of the summations is improved. We stress 
this is not just a solution to issues of limited accuracy. It also 
offers an improvement in MESA by providing a single solution 
methodology, network integration, that avoids the challenges 
of stitching together different solution methods. 

5.3. X-ray Burst Models and Adaptive Nets 

The new capabilities described above allow MESAstar to 
use large in-situ reaction networks (i.e. fully coupled to the 
stellar evolution rather than uncoupled co-processing). A 
demonstration is Type 1 X-ray bursts, a class of objects with 
unstable nuclear burning on the surface of a neutron star (NS). 
These burs ts are sensitive functio ns of accretion rate (Chen| 
et al.|1997|l, accretion composition (|Galloway et al.|2006|), the 
spatial distribution of burning on the surface of the NS (Bild- 


ample of such an equilibrium calculation is NSE, where a 
root-find for the neutron and proton chemical potentials is per¬ 
formed. Once these two chemical potentials are known, all the 
abun dances can be determined from nuclear Saha equations 
(e.g., IClifford & Taylerll 19651 IHartmann et al. 111985 IMeyer 

et^TJrP^ Seitenzam^^!pUn^ 


evolve the abundances and nuclear energy generation rate and 

sten 

1995 1 , the type of burning that occurs between bursts 

to replace it with equilibrium solution techniques. An ex- 

(Ga 

loway et al. 1200811 as well as possibly other conditions, for 


Odrzy wolek| I 

selves are efhi 


( jCumming & Bildsten|200l) l. Here we focus on a simplified 
model of constant accretion rate, where the burning occur s 
over the whole surface of the NS. GS 1 826-24 ([Tanaka T989 |), 


also known as the “clocked burster” (jUbertini et al.| 


19991, 


20121. Equilibriurn solution methods by them- 


rcient, robust, and inexpensive. 

However, combining reaction networks and equilibrium so¬ 
lution methods creates its own numerical issues, especially 
when the temperature and density are spatial and time depen¬ 
dent. Eor example, the temperature of a cell may start rel¬ 
atively low, move into quasi-static equilibrium (QSE) range 


provides an example of such a system due to its regular Type 
1 X-ray bursts. 

As material is accreted at the surface of a NS it is com¬ 
pressed and heats the underlying material. The a ccreted hy- 
drogen (from a low mass main sequence (MS) star (|Chen et al.| 
1997|l) burns via the hot CNO cycle. However, with high 
enough accretion rates the hydrogen will be accreted faster 
than the hot CNO cycle, which is limited by the jS-decay 
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Table 5 

Recurrence times of X-ray bursts. 


Model 

Accretion rate (10 ^Moyr 

Composition 

Recurence time (hrs) 

GS 1826-24 
rp_53 

3.00 

2 % metals 

4.0750 ± 0.0003 

1.5 ±0.10 

rp_153 

3.00 

2 % metals 

3.3 ± 1.80 

rp_3Q5 

3.00 

2 % metals 

3.2 ± 0.07 

rp_305 

3.00 

2 % ‘'•n 

3.0 ±0.07 

rD_3Q5 

2.40 

2 % metals 

4.1 ±0.30 

Heger et al. 2007 

1.17 

2 % ‘^N 

5.4 ±0.10 
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Figure 21. Kippenhahn plot during two X-ray bursts for the rp_3S5 net with 
the solar metallicity accretion model. The x-axis values are times relative to 
the peak of each burst, note the non-linearity of the scale. The y-axis values 
are the column depth and the color coding shows the temperature of the NS 
envelope. The dashed contours show the extent of the convective regions. 



Figure 22. The folded burst profiles for the different nuclear networks as 
compared to GS 1826-24 for an accretion rate of 3 X 10“* Mq yr“’ with 2% 
metals with a solar composition. Three rp network models are shown and 
one of the adaptive net models. The insert shows a zoom in of the first 30 s 
during the burst. 


timescale (of order minutes), can process the material. The 
accreted helium ignites unstably in a hydrogen rich environ- 
ment, allowing rapid proto n (rp) captures onto seed nuclei 
(Wallace & Woosley|1981[). This process forms nucle i along 


the proton d rip line up to and beyond the iron group (Schatz 
et al. II 1999 i, peaking at *°^Te, when g-decays prevent hea~ 


ier elements from being formed (|Schatz et al.||200T| |Fisker| 


et al. 20081. Once the burst begins, convection will com¬ 


mence, mixing the freshly burnt material with the ashes of 
previous burning episodes ( [Weinberg et al.|2()06 1 . 

GS 1826-24 has been studied by t he Rossi X-Ray Timing 
Explo rer (RXTE) over several years ( [Galloway et al. 2004[ 
[2008 | l. The bursts showed a decrease in the recurrence time 
between bursts, from 4.1 hr in 2000 to 3.56 hr in 2002, though 
during each observational epoch the bursts were consistent 
with each other. Based on the ratio of the burst energy to the 
persistent flux, it is assumed that the bursts are powered by 
hydrogen burning of solar metallicity material. 

We model the NS envelope using inner boundary con ditions 
for mass an d radius of Me = 1.4 Mq and /?e = 11.2 km (|Heger| 


et al.|2007|l, implying a gravitational redshift of 1 -H z = 1.26. 
ITie base of the envelope is composed of an inert layer that 
does not undergo reactions. The lumino sity at the base of th e 
envelope is settoL = 1.6x10^"^ergs s“' (Woosley et al. 2004|. 
We base our nuclear networks on the 304 species rp. net net- 



Figure 23. Folded burst light curve for the rp_305 net, with a solar 
metallicity accretion composition, shown for M = 2.4 x 10“^Moyr“^ and 
M = 3x 10“^ Mq yr“^, normalized to the peak flux measured for GS 1826- 
24. The insert shows the first 30 s of the burst. 
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work of Fisker et al. ( |2008 1 , which includes proton rich iso¬ 
topes up to ‘‘"Te. Isotopes above ®®Zn, which is the peak 
isotope in the mesa_2Q4.net, are included due to the proton 
captu res possible on hig h-Z isotopes during the peak of the 
burst ( [Fisker et al.|2006l l. We also include the effects of rota¬ 
tional mixing by setting a minimum amount of mixing in the 
NS envelope. This mixing, while having a physical motiva¬ 
tion (jPiro & Bildsten|2007 Keek et al.|2009|l, is there primar¬ 
ily to improve the convergence of MESA models by smoothing 
out the compositional gradients that form in the ashes of pre¬ 
vious bursts. We include the post-Newtonian correction to 
correct the local gravity in each cell for GR effects. During 
the burst we allow the accretion to continue. 

Our results are compared to the RXTE observations of GS 
1826-24 over bursts 9-20. Time resolv ed spectra were binned 
during the bursts’ rise time and decay (Galloway et al.|2008' 


jZamfir et al.|2012] l. Data output by MESA is not GR tirne cor 
rected, thus we set the burst times to be f' = f(l H- z), and 
average multiple bursts to produce a scaled light-curve. 

Figure l2T] shows the temperature profile during two X-ray 
bursts, forme rp_3Q5 net, accreting solarmetallicity material 
at a rate of 3.0 x 10“® Mq yr“'. At f' ^ -10 s the envelope ig¬ 
nites material and drives the formation of the first convection 
zone. This zone expands outwards in the envelope mixing the 
ashes from the bur ning at the base of th e envelope outwards 
to lower pressures ([Weinberg et al.|2006|l. As the burst decays 
the convection zone recedes outwards and by f' =» 150 s the 
envelope returns to its pre-burst temperature profile. 

We test three reaction networks, rp_53, rp_153 and 
rp_3Q5, each a modified form of that in [Fisker et al.[ ( [2008| l. 
Tablej^shows that increasing the number of isotopes in the re¬ 
action network increases the recurrence time and that all (for 
M - 3.0 X 10“^ Mq yr“*) have recurrence times 1 hr less 
than that of GS 1826-24. 

Figurej^and its insert show the folded light curves for each 
of the three rp reaction networks plus the GS 1826-24 obser¬ 
vations. The rise time is sensitive to the net, with the largest 
net matching the observed slow rise. The observed decay pro¬ 
file is also best matched by rp_3Q5. Burst to burst variations 
of the models decrease with increasing net size and can be 
further reduced by increasing the temporal resolution of the 
models. However, increasing the size of the net reduces the 
variation without having to increase the temporal resolution 
and also highlights the impact of MESA capability to include 
large nuclear networks. 

To achieve a better match to the GS 1826-24 recurrence 
time (see Table|^, we reduce the accretion rate to M - 2.4 x 


10 ®Mo yr '. However, Figure 23 shows that the light curve 


comparisons are not as good as iot t he higher M. _ 

GS 1826-24 was also modeled by[Heger et al.[([2007[) with 
accretion of hydrogen, helium and 2% ^“^N. For comparison, 
we run a model with this same composition with the rp_3Q5 
net. Table shows that the recurrence time decreases slightly 
when accreting 2% ^'^N rather than 2% metals. The model 
with metal accretion is in better agreement with both the light 

curve rise and decay. _ 

We now explore adaptive nets ([Woosley et al.|2004[l, where 
we allow MESA to determine which isotopes (and reactions) 
are necessary by assessing the available reaction pathways 
for the most abundant isotopes. The network is constructed 
by first finding those isotopes with an abundance above a 
threshold, Xi;^eep. and then introducing those isotopes which 
are connected by adding or removing protons, neutrons, or a 


particles. That determination is made via the additional pa¬ 
rameters Xn (i.e. neutron reactions) and Xp (i.e. proton and 
a reactions) potentially re-adding isotopes removed with the 
initiating Xkeep threshold. 

Accreting solar composition material at M = 3.0 x 

10 ®Moyr ' we follow the model to the second burst, find¬ 
ing a recurrence time of 3.1 hrs, comparable to that from the 
rp_3Q5 net (Table [^. The adaptive net has a better rise time 
profile than the rp nets, while the rp_3Q5 net has a better fit 
to the decay. This gives us confidence that the rp_3Q5 net in¬ 
cludes all relevant isotopes which drive the X-ray burst and 
thus is a useful approximation. For suitable values for the 
sensitivity of the adaptive net, the net limits itself to =» 400 
isotopes between bursts, which increases to ^ 600 isotopes 
during the burst. Variations of a factor 100 in the threshold 
parameters only change the isotope count by at most 50 iso¬ 
topes and do not affect the final results. 

6. CORE-COLLAPSE SUPERNOVAE 

The capability of using large, in-situ reaction networks 
without the need for equilibrium or co-processing techniques 
was described in Sectionj^and applied to X-ray burst models. 
We extend our demonstration of this capability by first consid¬ 
ering pre-supernova models. We then combine the advanced 
burning development with the implicit treatment of shocks 
discussed in Sectionj^to core-collapse supernovae models. 



log (Pc/g cm") 


Figure 24. Evolution of Fc and pc in solar metallicity, non-rotating M,- = 15 
and 30 Mq pre-supemova models. The curves are calculated using an in-situ 
204 isotope reaction network. Locations of the core carhon, neon, oxygen, 
and silicon ignition are labeled, as is the scaling relation Tc oc pc, and the 
Eji/k^T K 4 electron degeneracy curve. Regions dominated by electron- 
positron pairs, photodisintegration, and rapid electron capture are shaded and 
labeled. 


6.1. Pre-Supernova Evolution without QSE or NSE 

Eigure 1^ shows the - pc evolution of M, = 15 and 
30 Mq models from the onset of carbon burning until iron- 
core collapse. These non-rotating, solar metallicity mod¬ 
els used the 204 isotope reaction network described in Sec¬ 
tion!^ and MESAstar’s “Dutch” mass loss prescription with 
77=0.K These models have a; 2200 cells on the main-sequence, 
=» 3500 cells as the star becomes a red supergiant, and =» 2300 
cells at the onset of core collapse. At core collapse the final 
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Figure 25. Radial velocity and Ye profiles at the onset of core collapse for 
the Mj = 15 Mq model. Dashed curves show the results using a 22 isotope 
network and solid curves show the results using a 204 isotope network. Both 
models are evolved from the pre main-sequence to the onset of core collapse 
with their respective reaction network. The vertical gray lines mark the mass 
of the iron core as defined by the Ye jump. 



m [Mq] 


Figure 27. Mass fraction profiles of the ten most abundant isotopes within 
the iron core at the onset of core collapse for the M/ = 15 M© model evolved 
with the 204 isotope network. The entire iron core is in NSE and the mass 
fractions adapt to the changing temperature, density and Fg- 



Figure 26. Thermodynamic profiles at the onset of core collapse for the 
M,=15 M© model. Dashed curves show the results using a 22 isotope network 
and solid cuiwes show the results using a 204 isotope network. Both models 
are evolved from the pre main-sequence to the onset of core collapse with 
their respective reaction network. The vertical gray lines mark the mass of 
the iron core as defined by the Ye jump. 


masses are Mf — 13.0 and 15.2 Mq. The curves fall below the 
Tc Pc scaling relation because the core becomes partially 
electron degenerate, as indicated by tracks crossing the Fermi 
energy Epf{ksT) ss4 curve. Evolution towards lower density 
at nearly constant temperature signals ignition of a nuclear 
fuel. 

Figure shows the radial velocity and Yg profiles at the 
onset of core collapse for the M, = 15 Mq model. Dashed 
curves show the results using a 22 isotope network and solid 
curves show the results using a 204 isotope network. Both 
models are evolved from the pre main-sequence to the onset 
of core collapse with their respective reaction network. The 
vertical gray lines mark the mass of the iron core as dehned 
by the Yg jump, which is m « 1.43 Mq for the 204 isotope 
model and m =» 1.59 Mq for the 22 isotope model. The infall 
speed has reached =» 1000 km s ' just inside these iron core 
locations. 


Figure|^shows the thermodynamic prohles at the onset of 
core collapse for the Mi-15 Mq model. Dashed curves again 
show the results using a 22 isotope network and solid curves 
show the results using a 204 isotope network. The vertical 
gray lines again mark the mass of the iron core as dehned by 
the Yg jump in Figure The impact of these differences 
remains to be explored. 

Figure shows the mass fraction prohles of the ten most 
abundant isotopes within the iron core at the onset of core 
collapse for the M, = 15 Mq model evolved with the 204 iso¬ 
tope network. Each isotope shown dominates the NSE com¬ 
position at some location within the iron core, although we 
stress that no NSE or QSE approximation was used; the same 
204 isotope reaction network was used throughout the entire 
model from the pre-main-sequence to the onset of core col¬ 
lapse. 

The most abundant isotopes in an NSE distribution gener¬ 
ally have an individual Yg that is within a small range of the 
local Yg. A small spread usually exists due to nuclear structure 
effects. For example, the dominant isotopes at the center in 
Figure 27 are "'^Sc and ''^Ca. These isotopes have individual 


Yg of dr?29 and 0.417, respectivel y; c ommensurate with the 
central Yg a; 0.428 shown in Figure 26] The dominant isotope 
changes as the NSE distribution adapts to the rapidly decreas¬ 
ing density prohle and increasing Yg prohle. All the isotopes 
in the iron core eventually become part of the compact rem¬ 
nant after the explosion. However, the thermodynamic and 
composition prohles near the mass cut depend on the prohles 
interior to the mass cut. 


6.2. Core-collapse Supernova Explosions 

The envelope shock tests described in Section |4~9] show that 
the hydrodynamic solver in MESA meets the basic require¬ 
ments for shock propagation in a star. The AGB star model 
was selected because of the well behaved conditions of its 
envelope — a density structure that is smooth and monotoni- 
cally declining towards the stellar surface, and a uniform com¬ 
position. 

Here we explore the more challenging conditions associ¬ 
ated with a strong shock born deep in the stellar interior of 
a massive star. We study the dynamics of such a supernova 
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Figure 28. Multi-epoch snapshots of the hydrodynamical simulation from energy injection mimicking core-collapse supernova. The initial model is a ISM© 
star at solar metallicity, evolved with mass loss but no rotation, and employing a nucleai* network of 22 isotopes. We show the density (top row), temperature 
(middle row), and velocity (bottom row), versus Lagrangian mass (left column) and radius (right column). In each panel, the solid line shows the MESA results 
and the dashed line the VID results. 


shock and the explosive nucleosynthesis that takes place in 
the wake of the shock during the first second. The yields from 
explosive nucleosynthesis depend on both the energy and the 
power (characteristic energy deposition timescale), while the 
dynamics of the shock are primarily dependent on the total 
energy deposited. 


The starting conditions for the explosion simulations are the 
two 15 Mq pre-supernova models discussed above; one for the 
approximate 22 isotope network and one for the 204 isotope 
network. 

6.2.1. Explosion Dynamics 
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Since the focus here is on the dynamics rather than nucle¬ 
osynthesis, we expedite the MESA simulation by using the 22 
isotope network for both the pre-supernova and the explosion 
phases. Before triggering the explosion, we remove the iron 
core of the red-supergiant star by placing the inner bound¬ 
ary of the grid at a Lagrangian mass of 1.75 Mg. We trig¬ 
ger the explosion by depositing 1.52x10^' erg at a constant 
rate during 1 s. The artificial viscosity is raised during the 
energy deposition phase (i.e., I 2 =0.01), when the shock is 
at small radii, and lowered in the subsequent evolution un¬ 
til shock breakout (i.e., I 2 - 0.003). Since the binding en¬ 
ergy of the envelope to be shocked is -3.2x10^“ erg at the 
time we trigger the explosion, this choice of energy deposi¬ 
tion yields a total energy at the end of the deposition phase 


of 1.2x10^* erg. This is generally considered a standard value 
for a core-collapse supernova. 

Figure shows that the development of the explosion is 
analogous to the tes ts us ing the (low-density) envelope of an 
AGB star in Section [4!9| but with significant quantitative dif¬ 
ferences. Here, the shock born at the edge of the iron core first 
travels through the dense CO-rich core. At the outer edge of 
the He-rich shell, the shock traverses a steep density gradient 
to enter the low-density H-rich envelope. Hence, the shock 
crosses regions with densities ranging from ~10^gcm“^ at 
the edge of the iron core down to 10 '°gcm“^ at the stellar 
surface. 

The radii of the innermost shells are initially very small 
since they lie at the outer edge of the iron core. Consequently, 


they suffer considerable cooling from expansion. Figure 28 
shows a drop in temperature from a few 10® K down to 
~ lO^K at ~ 1 day. In addition, the supernova shock splits 
into a reverse/forward shock structure when it encounters the 
density drop at the transition between the He-rich core and 
the H-rich envelope. The reverse shock is the new feature, 
absent in the envelope shock test, that causes a significant 
deceleration of He-core material. The conversion of kinetic 
energy into internal energy causes this inner material to heat 
up, erasing the cooling effect from expansion. The innermost 
layers, which travel the slowest, will be shocked last. These 
innermost zones can evolve to temperatures ~ 10"^ K. It is in 
these innermost regions at late times that the differences be¬ 
tween MESA and VID are the largest. The offset occurs in 
a region of relatively high density (~ 10“®gcm“^) and low 
temperature (~ 10^ K). The offset in temperature between the 
MESA and VID simulations at late times stems from a differ¬ 
ence in the equation of state for metal-rich regions. MESA 
accounts for ionization through the OPAL equation of state 
table for metallicities z< 0.04. For higher metal abundances 
where OPAL tables are unavailable, MESA currently assumes 
full ionization while VID solves for the ionization state of the 
gas. Note that such density/temperature regimes are normally 
not encountered in stellar interiors. For other quantities and/or 
locations/times, the agreement between MESA and VID is ex¬ 
cellent. 

We also note that in the MESA simulation, two small spikes 
appear in the temperature and density profiles at < 2.5 Mg at 
=» 1000 s after the energy deposition phase. This feature is ab¬ 
sent in VID because VID uses a much larger viscous spread 
when the shock is in the helium-rich core {R < Rq). One 
can reduce or eliminate such spikes by increasing the vis¬ 
cous spread, although this may visibly smear the shock when 
it crosses the H-rich envelope — the current choice seems a 
suitable compromise. 


In contrast to the envelope shock test, this supernova ex¬ 
plosion configuration raises the temperature by a factor of 
about ten. Consequently, because Prad/Pgas Ip, the 

post-shock material becomes completely radiation dominated 
(Prad ^ Pgas)- If we neglect the binding energy and the ki¬ 
netic energy of the post-shock material, the post-shock en¬ 
ergy is of the order of the explosion energy. We indeed 
find a good correspondence between the post-shock temper¬ 
ature computed by MESA and the temperature obtained from 
(Po/oF)*^^ (where a is the radiation constant, Pq is a fitting 
parameter, typically of the order of the explosion energy, and 
V - 47rR^^/3 is the volume within the shock radius Psj,)- As 
expected, we also find that the shock accelerates (decelerates) 
in regions where PshP^t, decreases (increases) outward. 


6.2.2. Explosive Nucleosynthesis 

Here we compare the shock nucleosynthesis results from 
the two independent codes, MESA and VID. The same ini¬ 
tial 204 isotope pre-supemova model was the starting point. 
Our first test case is a strong explosion triggered by inject¬ 
ing 1.57x10^^ erg for 0.05 s and within 0.02 Mg of the mass 
cut, which is positioned at the outer edge of the iron core at 
1.5 Mg. The exact choice of explosion energy, deposition time 
scale, and mass cut is not strictly relevant. 

Figures 1^ and 30 compare the mass fraction profiles of 


MESA with a 22 isotope network, MESA with a 204 isotope 
network, and VID with a 54 isotope network. The first com¬ 
parison at 0.0 s shows the impact of mapping from the pre¬ 
supernova 204 isotope network to the networks used in the 
shock nucleosynthesis test. The next comparison at 0.05 s is at 
the end of the energy deposition phase. The final comparison 
at 42.7 s is after explosive nucleosynthesis has completed. In 
all cases, the silicon-rich and oxygen-rich shells are strongly 
influenced by the explosion; the former primarily for the pro¬ 
duction of ^®Ni and the latter primarily for the production of 
^^Si and The ^®Ni yields at 42.7 s are 0.092 Mg for VID, 
0.087 Mg for MESA with 22 isotopes, and 0.096 Mg for MESA 
with 204 isotopes. 

Overall the agreement between MESA and VID on this strong 
explosion is very good. The small differences between MESA 
and VID in Figures|^and[^are due to the difference in map¬ 
ping procedures. 

MESA follows rules for mapping isotopes from one network 
to another network; If an isotope present in the old network is 
also present in the new network, then the abundance from the 
old network is copied to the abundance for the new network. 
Isotopes in the new network that are not present in the old 
network are initially given a mass fraction of zero. MESA then 
separately renormalizes classes of isotopes to have the same 
total mass fraction in the new network as in the old network. 
The classes are neutrons, hydrogen, helium, carbon, nitrogen, 
oxygen, and other metals. This procedure guarantees that the 
sum of the mass fractions in a given class will be the same in 
the new network as in the old network. VlD’s mapping proce¬ 
dure for isotopes is the following. For any isotope present in 
the VID network but absent in the MESA input the mass fraction 
is set to the solar metallicity value. When an isotope included 
in the MESA input is absent in the VID network, this isotope is 
left out in the VID simulation. After completing the mappings, 
the resulting composition is renormalized so that the sum of 
the mass fractions is unity. 

Our second test case is a lower power explosion. We inject 
1.326x10^' erg in 1.0 s over 0.05 Mg of the mass cut, which 
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Figure 29. Nucleosynthesis profiles of selected isotopes for the 
1.57x10^* erg energy deposition test case at 0.0s (top) and 0.05 s (bottom). 
The dashed lines show the MESA results with a 22 isotope network, the solid 
lines show the MESA results with a 204 isotope network, and the long dashed 
lines show the VID results with a 54 isotope network. 



Figure 30. Same as Figure [29] but at a time after all nucleosynthesis has 
completed. 


is also positioned at the outer edge of the iron core at 1.5 Mq. 
The total energy after the deposition phase is 10^' erg. 



Figure 31. Composition profiles of the eight most abundant isotopes at 1 s in 
the inner 0.4 Mq of the ejecta for MESA (solid) and VID (dashed) simulations 
of a 10^* erg explosion in the 15 Me model. The Si-rich and 0-rich shells 
have been influenced by explosive nucleosynthesis, the former primarily for 
the production of ^®Ni and the latter primarily for the production of ^*Si and 

32s. 


In Figure!^ we show the composition profiles for the eight 
most abundant isotopes in the inner 0.4 Mq at the end of 
the energy-deposition phase (i.e., at 1 s). The correspondence 
between MESA and VID is again very good. The ^®Ni mass 
fraction approaches unity - it would reach unity if we appre¬ 
ciably increased the power (see Figurej^for example). Some 
^^Ni is produced in the same region, while ^^Fe is synthesized 
in the layers immediately above. The ^®Ni yields at 1.0 s are 
0.0041 Mq for VID, and 0.011 Mq for MESA with 204 isotopes. 

This work shows that the power of the explosion has a a sig¬ 
nificant impact on the abundance profiles. In the high power 
explosion, the yield of ^®Ni is 10 times larger and the ^He is 
several orders of magnitude more abundant. The nucleosyn¬ 
thesis of the low power explosion is completed at end of depo¬ 
sition phase at 1.0 s, while nucleosynthesis in the high power 
case continues for =» 30 s. This sensitivity suggests poten¬ 
tially observable signatures between low and high power ex¬ 
plosions. In addition, the explosive nucleosynthesis that takes 
place in core-collapse supemovae is sensitive to the way the 
explosion is triggered. With the approach we use (fixed power 
during the energy deposition phase), we find that increasing 
the explosion energy (at a given power), the power (at a given 
explosion energy), or both alters the amount of mass burnt. 
Moving the mass cut deeper into denser layers considerably 
enhances the amount of burnt material but this material may 
fall back rather than be ejected. Moving the mass cut fur¬ 
ther out into lower-density regions may completely quench 
the production of ^®Ni, in favor, for example, of ^^Si. It is 
thus important to keep in mind that the piston or thermal ex¬ 
plosion trigger is artihcial and that the yields from explosive 
nucleosynthesis bear significant uncertainties. 

7. IMPROVED TREATMENT OE MASS ACCRETION 

Adding mass to a star requires a way to accurately and effi¬ 
ciently compute the thermal state of the freshly accreted ma¬ 
terial in the outermost layers. This is simplified by a hier¬ 
archy of timescales. Eor accretion at M, there are two im¬ 
portant timescales at a given location, m, the thermal time 
Tth - {M-m)CpT jL, where Cp is the specihc heat at constant 
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P, and the time to accrete this same layer, Tacc - {M - m)IM. 
Near the surface, L » CpTM, implying that rth ^ Tacc, so 
that these layers have ample time to relax to the thermal equi¬ 
librium configuration hxed by L (|Nomoto & Sugimoto|1977[ 


Nomoto|1982||Townsley & Bildsten|2004|l. 

cases where L arises solely from compression of mate¬ 
rial, such as a very rapidly accreting star of high M or an old, 
cold accreting WD, then L ~ CpT\,M, where T\, is the tem- 
perature at the degenerate/no ndegenerate transition in a WD 
(|Townsley & Bildsten||2004] l or of the core in a normal star. 
Even in these cases, the outer layers have T ^ Tj,, allowing 
the inequality Tth Tacc to hold. This also implies that the 
thermal state of the arriving material is unimportant, allowing 
us to safely use the approximation that material arrives with 
the same entropy as the photosphere, since material relaxes 
toward this on the very short Tth at the photospherej^ Even 
when M varies on short timescales, using an average accre¬ 
tion rate is a good approximation for computing the evolution 


of the interior layers due to their long Tth (Piro et al. 2005 
Townsley & Gansicke|20091l. 

The timescale hierarchy Tth ^ Tacc implies that the outer 
regions evolve nearl y homologously in the fractio nal mass co¬ 
ordinate q - rrijM |Sugimoto & Nomoto|1975'] l. Hence, the 
thermal profile (e.g. the run of T with P or p) of the outer 
layer is nearly constant in time even as fluid elements are com¬ 
pressed to higher pressures and have Tim) increase. More 
formally, Tiq) varies slowly in time near the surface, where 
(1 - ^) « 1. This motivates reformulating the Lagrangian 
based form of 


tgrav 


Ds (ds 

= -T — = -Ti¬ 

nt \dt 


(60) 


that is needed in the energy equation, dLldm - egrav + 


( |Sugimoto]|1970[ |Sugimoto & Nomoto||1975t Paper I), to 


version in the coordinate q. 


*-grav 


= -T 



j dinMjt) 

dt 


ds \ 
d\nq), 


( 61 ) 


Sugimoto & Nomoto (19751, and later works based on it, de- 
note the second term on the right hand side the “homologous” 
term. Physically it is the local loss of entropy in the fluid el¬ 
ement as it is compressed to higher pressure. They label the 
hrst term on the right hand side the “non-homologous” term. 
It arises from the much slower departure of the outer layers 
from simple homologous evolution on a timescale M/M. 

MESAstar includes the ability to have an inner inert core of 
mass Me- In this situation q - (m- Me)KM - Me) rather than 
the more typical m/M. For simplicity here, we will use Me-Q. 
Approximate homology holds in either case for 1 - o' 1. 


In P aper II, following the work of [Townsley & Bildsten| 
(|2004[), only the homologous term in Equation (|61|l was in¬ 
cluded in egrav in and near regions of newly accreted material. 
We also described the huge advantage of such an implemen¬ 
tation, as it allows the mass added per timestep to be much 
larger than the smallest cell mass near the surface, while main¬ 
taining accurate thermal profiles at low pressures. However, 
leaving out the non-homologous term can create a disconti¬ 
nuity in egrav at the location where the standard Lagrangian 
derivative. Equation (|60|l, begins to be used. We now de- 


A possible exception to this case is rapidly accreting pre-main sequence 
stars where the accretion shock is so o ptically thick that the m aterial’s entropy 
remains high and is advected inward jPalla & Stahler|1990^. 


scribe the improvement we have made to MESAstar so that 
it now includes both the homologous and non-homologous 
terms. Hence, the two forms of DsjDt are physically equiva¬ 
lent and there is no longer any discontinuity. 


7.1. Lagrangian and Homologous Regions 

The independent coordinates used for writing the time- 
dependent structure of the star are m and f, and for fluid ele¬ 
ments deep within the star at both timesteps, the conventional 
form of Equation ( |60l l is adequate. One numeric subtlety of 
accretion is that the derivative at constant m cannot be evalu¬ 
ated for material that is not present in the star at the beginning 
of the timestep. However, the simplification available when 
Tth ^ Tacc, manifest in Equation enables egrav to be eval¬ 
uated in the outer regions. 

When T and p are used as independent variables, we write 


cgrav 


= -CpT 


(1 - '^hAXt) 


d\nT 

dt 




(62) 

in the interior of the star, as in Paper I, and near the surface 
we choose to write, using Equation (|6T]i, 


't^grav “ ‘“grav, 


^srav.nh ^srav.h , 


(63) 


where 

^grav.nh “ CpT 

and 


/(91n7’\ /51np\ 

(l-VadTT)(^^j -VadA-pJ^j 


CpTGmM 

egrav,h - ' 


(64) 


(65) 


Here Vp = c/lnT/c/lnP is the T-P profile in the star, 
and we have used the thermodynamic derivatives Vad = 
id\nTld\nP)s, Xt - {dPldT)p, and;^p = (dPldp)T- There 
is also a transition region where a weighted combination of 
these forms is used, with weights varying linearly in m. 


Lagrangian transition homologous 


constant m 

modified m 






\\///K^ 


\^6M^ 


Mass coordinate m m = M 


Lagrangian transition homologous 


modified q 

constant q 

\ 














center Homology coordinate q q=i surface 

Figure 32. Illustration of the MESAstar mesh in both m and q for a timestep 
in which mass is added to the star over the time interval to / 2 - Vertical 
and slanted lines indicate cell boundaries. Cell size is exaggerated; there can 
be many cells in the newly added material. Three regions are chosen in the 
process of expanding the mesh for the added material, an inner Lagrangian 
region, an outer homologous region, and a transition region. 
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Placement of the transition is related to the mesh. The 
MESAstar mesh structure is unchanged from that discussed 
in Paper I and Paper II. An illustration of the mesh regions 
and an indication of the behavior of cell boundaries during 
accretion is shown in Figure [32] In the case of mass loss, 
analogous operations are pertomed; here we focus only on 
mass gain. 

A mass 6M - M{t 2 - fi) is added to the star from time fi 
to f 2 . Sizes are exaggerated for clarity; there are generally 
many zones in each region, possibly hundreds in the newly 
added material. The diagram is shown in both mass coordi¬ 
nate, m, and homology coordinate, q. Before each timestep, 
MESAstar adjusts the initial mesh resolution by splitting or 
merging cells based on local gradient conditions, producing 
an adjusted resolution mesh, which we show at t\. No sim¬ 
ulation time elapses during that process. When constructing 
the mesh that will represent the star at f 2 , MESAstar divides 
the new mesh into three regions: an inner Lagrangian region, 
in which the m boundaries of each cell are preserved during 
the timestep, a transition region, and an outer, homologous 
region, in which the q boundaries of each cell are the same 
across the timestep. The result of this operation is an ex¬ 
panded mass domain shown at f 2 in Figure 

Time derivatives appearing in the equations for physical 
evolution are then estimated using first order differences. In 
the Lagrangian mesh region, the finite difference form of 
(dfdt),„ involves a simple same-cell difference. Similarly, 
in the homologous mesh region, the finite difference form 
of {dldt)g involves a same-cell difference. In most cases, 
by design, these same-cell differences are for values whose 
changes, e.g. ^InT, are directly available from the iterative 
solution of the new structure, allowing us to avoid the nu¬ 
merical problems inherent in subtracting two almost identical 
numbers. In the transition region both m and q coordinates 
of cells have been modified, so we cannot do a same-cell dif¬ 
ference for either {dldt)„t or (dldt)^. Instead, we interpolate 
values from the model at the start of the step to corresponding 
locations in m or q at the end of the step. 

A smooth and accurate value for egrav in the transition re¬ 
gion is important. To ensure this, the location of the transition 
region is selected to reduce the differences between the con¬ 
stant m and constant q forms of the time derivative and main¬ 
tain accurate finite differencing. As a simple mechanism to 
control these, we limit, in units of cell size, the offset in the 
interpolation used to tran slate locations from the beginning to 
the end of the timestep ( [Miles et al.||2015] l. Using the cell 
size implicitly takes advantage of the limits imposed by mesh 
controls on the maximum possible magnitude of cell-to-cell 
changes in key variables, including the variables of interest 
for egrav time derivatives. 


7.2. Testing 

In order to demonstrate that the interface between the outer 
homologous region and the inner Lagrangian region provides 
a smooth profile that is independent of timestep size, we have 
repeated the test shown in Paper II section 5.3. This test in¬ 
volves accretion of solar composition material onto a WD 
at IO^'^Mq yr“'. We use the same starting model as in Pa¬ 
per II, which was produced by accreting hydrogen-rich ma¬ 
terial through several hydrogen shell flashes on a Q.6Mq WD 
with an initial core temperature of about 10^ K. As accre¬ 
tion proceeds, the total accumulated accreted mass, Macc, in¬ 
creases up to a maximum which causes the hydrogen flash and 


nova runaway, Mign. 


log(P/ergcm (approximate) 



log(l - q) 

Figure 33. Profiles spanning the Lagrangian-homologous grid transition re¬ 
gion in a 0.6Mq WD with a core temperature of 10^ K undergoing accretion 
of solar composition material at a rate of 10“^^Mo yr“^. At the time shown 
Macc/^ign = 0-2. Shown for comparison are the results of the treatment dis¬ 
cussed here for a 5000 yr timestep and a 1000 yr timestep and the treatment 
discussed in Paper 11. For the current treatment the transition from homolo¬ 
gous to Lagrangian is indicated by two open circles at either end of the region. 
The location 6M in from the surface, the base of material newly accreted this 
step, is indicated by the triangle. 


Profiles near the transitions region are shown in Figure 
at the time when Macc/Afign = 0.2, which had the most se¬ 
vere discontinuity in Paper II. The first panel shows and 
the second panel shows Anr^pHp x egrav, where Hp is the 
pressure scale height. This is the amount of energy being re¬ 
leased due to the T DsjDt term in the energy equation within 
a scale height, and has units of luminosity. We have chosen a 
timestep of 5,000 yr, which places the homology-Lagrangian 
transition region in a similar place to the location of the dis¬ 
continuity in Paper IT 

The orange curve in Figure [^ was computed using MESA 
version r4664, as used in PapetTl. This displays the disconti¬ 
nuity due to using egrav from Equation ( [62| i and Equation ( |6^ 
with only the homologous term and a transition point a fac¬ 
tor of 5 deeper in pressure than 6M. The black curve shows 
the same simulation with the same timestep for the treatment 
discussed here, in which egrav,nh is included and the transition 
region, indicated with solid dots at either end, is placed as de¬ 
scribed in Section [7T| We see that, away from the transitions, 
egrav is unchanged from the values found in Paper II, in which 
only the homologous term was used in the exterior. We also 
show the result for a timestep of 1,000 yr is indistinguishable 
on this plot; the egrav profile differs by a fraction of a percent 
at the edges of the transition region, and less elsewhere. In 
the current treatment the profiles are independent of timestep 
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size. 


8. WEAK REACTIONS 


An erratum for this paper was published as Paxton et al. 
(20161. The errors reported there have been coirected in this 


version of the paper. 

The rates module provides weak reaction rates for hun¬ 
dreds of isotopes. By default, when atoms are fully ionized, 
these rate s are based (in order of precedence) o n the tabu- 
lations of |Langanke & Martmez-Pinedo| (|2000|), |Oda et al. 
(|1994|), an^Fuller et al.f^l985f^ ITiese tables span a wide 
range of density and temperature, 1 < logOoTe/g cm 3) < 11 
and 7 < log(7’/K) < 10.5, but are relatively coarse, with 11 
points in the pY^ dimension (AlogpTe = 1) and 12 points in 
the T dimension (A log T =» 0.25). 

These grids include the thermodynamic conditions where 
the electrons are degenerate and relativistic, which are real¬ 
ized for example in massive white dwarfs and cores of in¬ 
termediate mass stars. Under these conditions, the rates of 
electron-capture and beta-decay reactions are sufficiently sen¬ 
sitive to density and temperature that they can change by tens 
of orders of magnitude between adjacent points in these ta¬ 
bles. Linear or cubic interpolation cannot accurately repro¬ 
duce the value of the rate between the tabulated points. 

The difficulty of interpolating i n coa rse rate tabulations 
was discussed by [Fuller et al.| ( |1985| l, who proposed a 
physically-motivated interpolation scherne, hereafter referred 
to as Fuller, Fowler, and Newman (FFN) interpolation. Their 
procedure assumes the rate has the form given by a single 
transition between the parent and daughter nuclear ground 
states. However, the true rate may be dominated by allowed 
transitions to or from excited states in the parent or daugh¬ 
ter nucleus. This is almost always the case when the ground 
state to ground state transition is highly forbidden. The spe¬ 
cific transition that dominates the rate may change over the 
range of thermodynamic conditions covered by the table. The 
FFN interpolation method does not account for these compli¬ 
cations. 

Figure compares the results of the interpolation meth¬ 
ods descrmed in the preceding paragraphs with the on-the- 
fly approach that we have implemented in MESA and will be 
described here. It shows the electron-capture rate on ^^Mg 
and beta-decay rate of ^"^Na at fixed temperature. Linear in¬ 
terpolation of these coarse tables fails to reproduce the rapid 
variation in the rate. The FFN interpolation method produces 
curves with characteristic shapes more similar to the true rate, 
but because the Q-value is that of the ground state to ground 
state transition and not that of the transition that dominates 
the rate, the density dependence is not correct in detail. 

In recent years, a number of authors have discussed the 
importance of well-sampled weak rates in capturing the in¬ 
fluence of these processes on stellar evolution. This can be 
achieved by generating denser tables for the specific reactions 
of in terest or by using analytic approximations to the rates 
(e.g., Toki et al.|2013| Martmez-Pinedo et al.|2014| ). We now 
present a capability by which MESA can calculate weak reac¬ 
tion rates on-the-fly from input nuclear data. This removes 
the potential for interpolation artifacts. It also enables easy 
experimentation in cases where the input nuclear data may 
not be well-measured. We begin with an overview of how we 
calculate these weak rates and illustrate their utility and a few 
applications. 

8.1. Calculation of Weak Rates 


o 



to 


Figure 34. The top panel (bottom panel) shows the rate of electron capture 
on (beta decay o f ^‘^Na) as a functi on of density at a fixed temperature 
of log(7’/K) = 8.6. The|Oda et al.H1994) tabulated points are shown as black 
dots. The dotted line shows the result of using linear interpolation between 
the tabulated points. The dashed line shows th e result of using the physically- 
motivated interpolation method suggested by|FuIIeretak|(T985). The solid 
line shows the rate calculated using the on-the-tly rate calculation capability 
of MESA documented in this section. Slight differences between points and 
the line are due to differences in the input nuclear data. 
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Consider two nuclei A = (Z,N) and B = (Z - l,N + 1) that 
have two states connected by an electron-capture transition 

A H- B -H Vg (66) 

and beta-decay transition 

B —> A + e + Vg . (67) 


The energy difference between the ground states can be writ¬ 
ten as 


_((Ma-Mb)c^ 
\(Mb - Ma)c^ 


for electron capture, 
for beta decay. 


( 68 ) 


where and Mg are the nuclear rest masses of the ground 
states. The total energy difference between any two states can 
be written as 


Qij^Qs + E,-Ef, (69) 


where £, and Fy are the energies of the initial and final states 
measured relative to the ground state. For the transitions that 
we consider here, Qg <0 and Qij < 0 for electron capture and 
Qg> 0 and Qjj > 0 for beta decay. 

In this section, we use J to represent the nuclear spin. We 
work in the allowed approximation, which neglects all total 
lepton angular momentum (L = 0). This restricts us to Fermi 
transitions, where the total lepton spin is 5 =0, and therefore 
the initial and final nuclear spins are equal (7, = Jf), and 
Gamow-Teller transitions, where 5 = 1, and therefore 7, = 
Jf, Jf ± I (excluding 7, = 7/ = 0). In both cas es, there is no 
parity change: n/nf = H-1 (e.g., Commins|1973|l. 

The total rate of the process (electron capture or beta decay) 
is the sum of the individual transition rates from the f-th parent 



































33 


state to the j-th daughter state, Aij, weighted by the occupation 
probability of the /-th parent state, p, . 

(70) 


where (ft) is the comparative half-life and can be either mea¬ 
sured experimentally or calculated from theoretical weak- 
interaction nuclear matrix elements. The /-sum is over all 
parent states and the y-sum is over all daughter states. The 
occupation probability is 

^ (27, + l)exp(-^£,) 

2^(27, + 1) exp(-^£,) ’ ^ ^ 

where we define (3 - (kfT)~^. The quantity O is a phase space 
factor which depends on the electron chemical potential pe 
(including the electron rest mass), on the temperature T, and 
the energy difference Qij. The value of O for electron capture 
is 


_ eXp(TOZ) r“ £’e('£’e + Qij) 

(mec2)5 J_Q,^ 1 H- exp[/3(£'e - jUe)] 



Figure 35. The electron capture (solid lines) and beta decay (dashed lines) 
rates of the ^^Na-^^Ne Urea pair as calculated by MESA, using the on-the-fly 
methods described in this section. The value of log(7/K) is shown next to 
each electron capture line; the beta decay line of matching color is at the same 
temperature. The rates vary rapidly, with both temperature and density, near 
the threshold density, which is roughly in the center of the plot. 


where a is the fine structure constant. For beta decay it is 


O/? - 


exp(7rQrZ) 


rQij 


El(E, - Quf 


-dEe 


(mec2)5 1 -H exp[-/3(£'e - jUe)] ' 

Similarly, the total rate of energy loss via neutrinos is 
meAln2„,^ ^ , 

(ftfj 

The value of'T for electron capture is 


(73) 


(74) 




exp(;7raZ) 


f 


(m^c^)^ j^Q.. 

and for beta decay it is 
exp(TOZ) 


EjiE, + Qij)^ 

1 H- exp|;8(£'e - Ee)] 


dEe 




(nif^c^)^ 


r'da 


El(E, - Qijf 

+ exp[-y8(£'e - Ee)] 


dEe 


(76) 


In order to implement these equations in MESA, we rewrite 
the integrals in terms of Fermi-Dira c integrals, following ap¬ 
pendix A of ISchwab et al.| (|2015|l. MESA implements fast 
quadrature routines to evaluate integrals of this form. Each 
time a weak rate is needed, it is calculated on-the-fly. We dis- 


cuss the computational cost of this procedure in Section 8.3 
Assuming thermal equilibrium, the energy released by 
weak reactions depends only on total reaction rate, total neu¬ 
trino loss rate, energy difference between the nuclei, and the 
electron chemical potential. Therefore, the total specific heat- 


of electron capture and beta decay and the corresponding en¬ 
ergy generation rates. Typically only a few low-lying states 
and the transitions between them are needed. As an example. 
Figure [^shows the rates for the ^^Na-^^Ne Urea pair. 

8.1.1. Coulomb Corrections 

In a dense plasma, the electrostatic interactions of the ions 
and electrons introduce corrections to the weak rates relative 


to those which assume a Fermi gas of electrons and an ideal 
gas of ions. Our treatment of these effects, which is pre- 

sented in appendix B of Schwab et a 

.|(2015|l, is similar to 

(/6) appendix A of|Juodagalvis et al. 

PJTU 

ecay change the ion 

Since electron capture and beta c 


charge, the Coulomb interaction energy changes the energy 
difference between the parent and daughter nuclear states. To 
calculate this shift, we use the e xcess ion chemical potential 
from Potekhin et al.|(2009|l. We incorporate this effect 




by shifting the value of Qij, as defined in Equation (|69]l, by an 


amount AE = 


This shift, Q'.. - Qij + AE, 


-^ex,parent /^ex,daughter ■ 

enters the calculation of the phase space factors and the en¬ 
ergy generation rates. 

The electron density relevant to the reaction rate is not the 
average electron density, but rather the electron density at the 
position of the nucleus. This correction is accounted for as a 
shift in the value of the electron chemical potential that enters 
the phase spa ce factor, /i' = + Vj- Values of Vj have been 

calculated by Itoh et al. (20021. This correction does not en¬ 


ter the energy generation rates because it has not changed the 


ing rate from a reaction is 


energy cost to add or remove an electron. 


\iQg A^e)'^ec ^v,ecj » 

(77) 

8.2. Applications 

nB 

[(Qg - \ , 

(78) 

When pe < 121, only the few electrons in the tail of the 
Fermi-Dirac distribution have sufficient energy to overcome 


where and ng are the number densities of the species un¬ 
dergoing electron capture and beta decay, respectively, and p 
is the total mass density. 

Therefore, given a list of nuclear levels and the (/f)-values 
for the transitions between them, MESA can calculate the rates 


the energy gap and capture on A to form B. Thus, the rate of 
electron capture is small compared to the rate of beta decay, 
and so isotope B is favored in the equilibrium. When > 
121, there are only a few unoccupied states available to accept 
the energetic electron from the beta decay. This final state 
blocking means the rate of beta decay is small compared to 
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the rate of electron capture, and so isotope A will be favored 
in the equilibrium. 

The shift in this equilibrium can have profound conse¬ 
quences when it occurs in stellar interiors. It modifies the 
composition, reduces the electron fraction, and alters the ther¬ 
mal state of the plasma. We now discuss two applications of 
our on-the-fly treatment of the weak rates: the Urea process 
and accretion-induced collapse. 

8.2.1. Urca Process 

When the ground state to ground state transition is allowed 
(odd nuclei), the rates of electron capture and beta decay are 
both significant when jUe a; \Q^\. Since each reaction produces 
a neutrino which free-streams out of the star, this can lead to 
significant cooling. With a total number density of an Urea 
species tiu - ha + ns, assuming the abundances are given by 
the detailed balance condition riA^-ec + nB^p — 0, the volumet¬ 
ric neutrino cooling rate will be «t/C, where 


C = 


^v,ec + EyfiA ei 


(79) 


Ap + Tec 

In the limit k^T I2I, the maximum value of the Urea cool- 


ing rate at a given temperature has a simple form (e.g., Tsuruta 
|& Cameron|1970| l 


C - 

^ max 


77rMn2 


60 \{ft)p + {ft\ 



4 

( Q ) 


1 1 

[m^c^ I 


exp(7raZ) . 
(80) 


Well-sampled rates such as those shown in Figure 35 are nec¬ 
essary to reproduce the correct Urea cooling rates^We illus¬ 
trate this in Figure which shows Umax for the ^^Na-^^Ne 
Urea process for temperatures 10® - 10® K. The circles show 
the results using the on-the-fly treatment described in this pa- 
per; the squares show the results using the coarse tables of 
Oda et al. ( | 1994 1. The dashed line shows the cooling rate 
expected from Equation ( |80| ) which is in excellent agreement 
with the results of the on-the-fly method. The Urea cooling 
rates calculated from interpolating in coarse tables severely 
underestimate the true cooling rate when k^T <si \Q\. 

Thus, when the Urea process is important, well-resolved 
weak rates are necessa ry to correctly cap t ure the temperatur e 
evolution of the co re ( |Toki et al.||2013[ [Jones et al.||2013| l. 
[Jones et al.j ( [201311 used M ESA r3709 along with a denser 
table describ ed in [Toki et al^ ( [2013[ l to do their work. The [Toki[ 
et al.|(20T3]l table is not publicly available, so to reproduce the 


results of [Jones et ar] ( [2014[ ) we save a model of an 8.8 Mq star 
at log(pc/gcm~^ ) - 8.95 from ou r run with MESA (version 
r3709) using the Jones et ar] (2014l inlists. We then load this 
model into a newer MESA version (r7503) that has access to the 
on-the-fly weak rates and evolve this model using a network 
with only the Urea process reactions. During this phase other 
nuclear reactions are not important to the central evolution. 

Figure]^ shows the central temperature and density of the 
core. The solid lines show the evolution using the on-the-fly 
rates, the dashed lines show the results when interpolating in 
coarse tables. The drops in temperature at log(pc/gcm“®) 

9.1 and « 9.25 correspond to cooling from the ^^Mg-^^Na and 
^®Na-^®Ne Urea pairs, respectively. The corresponding shifts 
in composition can be clearly seen in the lower panel. These 
results demonstrate the importance of densely-sampled weak 
rates to the evolution of the core. 


8.2.2. Accretion-Induced Collapse 
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Figure 36. The effect of the interpolation method on the Ure a pro cess cool¬ 
ing rates. The circles show the maximum value of C (Equation |79^ calculated 
using the on-the-fly methods disc ussed in thi s sectio n; the squares show the 
results using the coarse tables of|Oda et al.Hl994). Interpolation in these 
coarse tables severely underestimates the Urea cooling rates at low tempera¬ 
tures. Th e das hed line shows the expected value of the cooling rate given by 
Equation {SOj . 
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Figure 37. The top panel shows the evolution of and pc in an 8.8 Mq 
star. The bottom panel shows the central and ^^Na mass fractions. 

The solid lines show the evolution using the on-the-fly rates, the dashed lines 
show the results when interpolating in coarse tables. The locations of the 
changes in mass fraction match the locations of cooling in the top panel. This 
demonstrates the importance of densely-sampled weak rates to the evolution 
of the core. 


When the ground state to ground state transition is forbid¬ 
den (even nuclei), the first transition to become significant is 
typically an allowed transition into an excited state. In these 
cases, the beta decays from the daughter ground state are 
blocked and decays from daughter excited states are strongly 
suppressed by the Boltzmann factor. Therefore significant 
cooling via the Urea process does not occur. Instead, since 
the captures are preferentially to an excited state, significant 
heating occurs via gamma-ray emission as level populations 
relax to a thermal distribution. 

Two important capture chains occur in oxygen-neon- 
magnesium (ONeMg) cores: ^^Mg ^ ^'^Na —> ^'^Ne and 
^^Ne —» ^®F —» For these sequences of captures, the 

excess electron energy is thermalized. These are the key 
reactions in electron-capture supernovae and the accretion- 
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Figure 38. The evolution of a cold ONeMg WD towards AIC. The top panel 
shows the evolution of pc and Tc. The bottom shows the central and 

^®Ne mass fractions. The solid lines show the evolution using the on-the-fiy 
rates; the dashed lines show the results when interpolating in coarse tables. 


induced col lapse (AIC) of ONeMg white dwarfs (e.g., Miyaji 
|et al. [T980| ). As the degenerate core approaches the Chan- 
drasekhar mass, the electron captures remove the pressure 
support and heat the plasma. Figure 1^ shows the evolution 
of a cold ONeMg WD {Xq = 0.5, Xnc = 0.45, Xug = 0.05) 
accreting at M = 10“^ Mq yr“'. The solid lines show the evo¬ 
lution using the on-the-fly rates described in this section; the 
dashed lines show the results when interpolating in coarse ta¬ 
bles. When using the coarse tables, the electron captures on 
^'^Mg do not occur until approximately a factor of two larger 
density. At this greater density, the energy deposition from 
each capture is higher and this leads to a large temperature 
change due to the A = 24 captures alone. In contrast, the 
on-the-fly rates show the behavior demonstrated in previous 
studies of this evolution that did not use sparse tables (e.g., 
Miyaji & Nomoto|1987|l: the A = 24 captures heat the plasma 
and accelerate the contraction; the A = 20 captures, due to the 
higher ^‘'Ne abundance and a higher energy release per cap¬ 
ture, cause a thermal runaway and the formation of an oxygen 
deflagration (Schwab et al.|2015]l. 


8.3. Guidelines 

MESA provides the nuclear data used in the ca lculation of 
the reactions specifically discussed in this section (Tilley et al. 


1998[ |Firestone|2007a|b[ |20091 IShamsuzzoha Basunia|20lTf 


Martlnez-Pinedo et al. 2014| l. To consider additional reac- 
tions, a list of nuclear levels and ( /f)-values must be specified. 

The expressions in Section [8T| assume degenerate, relativis¬ 
tic electrons. As /Zg increases, additional transitions to higher 
energy states of the daughter nuclei and must be included. At 
higher temperatures, excited states of the parent nucleus will 
begin to be thermally populated and captures or decays from 
those states and must be included. At temperatures and den¬ 
sities where the composition approaches NSE, these methods 
are particularly inappr opriate, as it is necessai y to consider 
large pools of isotopes ( |Juodagalvis et al.|2010| l. 

9. CHEMICAL DIEEUSION 

mesa’s early implementation of micr oscopic element dif fu- 
sion incorporated the approach used by|Thoul et al.|(]1994|l in 
their seminal work on understanding the sedimentation of he- 
lium in the solar interior. The fundamental starting point for 


this treatment of diffusion is the Boltzmann equation with the 
assumption of binary collisions where the particle’s mean free 
path is much larger than the average particle spacing. Thi s 
formalism, encoded in the Burgers equations ( |Burgers|1969| ), 
assumes that ions interact with an effective potential that gov- 
erns isolated interactions between only two particles at a time. 
Eor more strongly coupled plasmas, as F ^ e^HAionk^T) ex¬ 
ceeds unity (where dion = (3/47rnion)^^^ is the mean inter-ion 
spacing, and nion is the total ion number density), it is no 
longer clear that this assumption remains valid. L ater updates 
to MESA incorporated the work of|Hu et al.|(|201 1)1 on radiative 
levitatio n and incorpora t ed the resistance coefficients calcu¬ 
lated by [Paquette et al.| (|1986|) for approaches to the denser 
plasma regime as F —» 1. 

Here we describe MESA’s current implementation of chem¬ 
ical diffusion and then discuss the path forward for diffusion 
implementations in the F > 1 regime, needed for accurate 
studies of diffusion in the interiors of white dwarfs or surfaces 
of neutron st ars. Recent theoretical work in this strongly cou¬ 
pled regime (Baalrud & Daligault||2013 2014[ Beznogov &| 


Yakovlev 2014p provides support for a future update of MESA. 


9.1. Current Methods in NESAstar 

We now describe the formalism and assumptions underly¬ 
ing the approach to diffusion currently present in MESA. This 
is followed by a discussion of the framework for numerical 


implementation of this formalism provided by Thoul et al. 


(19941 and key modifications present in the current version 
of the MESA diffusion routine. 

9.1.1. Burgers Equations and the Low Density Limit 

The Burgers equations for diffusion in an ionized plasma 
are derived using the Boltzmann equation for the distribution 
function F'j(x, t) for particles of type s 


dt 


V—1 dF s 


dxi 


■y fsi dFs 


dt 


(81) 


collision 


where Xi are the components of the position vector, are the 
components of the velocity vector, fsi are components of the 
forces on particles of type s, and is the mass for those par¬ 
ticles. Throughout this section, the indices s and t refer to 
particle species, while i and j are used to index other quanti¬ 
ties such as spatial components of vectors. _ 

Burg ers adopts the 13-moment approximation due to Grad| 
( |1949[ ) as a closure scheme for taking moments of the Boltz- 
mann equation. Burgers also assumes an approximately 
Maxwellian distribution function 


Fs = 


r3/2«3 


exp 


(1 + fs). 


(82) 


where Qs - (2kB7’/mj)'^^, Csi - - Usi represents the com¬ 

ponents of the deviation of the velocity from the mean flow 
velocity Uj of the species, and 












(83) 


is the small deviation (fs ^ 1) from the Maxwellian distribu¬ 
tion. The coefficients Bsij and Csi are defined such that the 
distribution function has a total of 13 free parameters corre- 
sponding to the 13 moments of the closure scheme (see|Burg- 
|ers| 1969 1. 
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Burgers derives the collision integrals (5®) and cross- 
sections (2®^) that result from taking moments of the right 


hand side of the Boltzmann equation 

/^CC 

: I (1 - cos'Xst)bdb, 

Jo 


5® = 27r 


= I «?vexp(-^ 


47r 


rV2 


f 


St / Ct 


1 »2 j+3 
^ c(0 

2 j+A SI ’ 


(84) 


(85) 


where a^, - Ik^TIfist, l^st — msintHms -i- m,), v represents 
the relative velocity of colliding particles, and the angle of 
deviation Xst is a function of both v and the impact parameter 
b that depends on the physics of the two-particle interaction 
between colliding particles in the gas. Burgers then defines 
the dimensionless coefficients Zst, z'^,, z",, and z"', along with 
resistance coefficients (Ks,) in terms of the collision integrals: 




( 11 ) 


= -d 

y(13) ,y(ll) _ ^ _ 
^st / ^st ~ ^ 

y(22) /y(l 1) _ ff 
^St l^st ~ ^SV 
(11) _ /// 
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- Zst), 
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-zrZst ■ 


( 86 ) 


2^ St, 


y(23) /yo i; _ t 

2^ St 1St ~ s-s 


In the “single-fluid picture” the diffusion velocities are de¬ 
fined with reference to the mean velocity of the gas as a whole 
(u), rather than with respect to the mean species velocity (Uj): 


=-r 

ns J 


- d^^tFs, 


m^-Yps»s 


y/s - Uj - u. (87) 


Burgers defines residual heat flow vectors 


2nsksT 
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• u| V,, 
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( 88 ) 


As shown in section 18 of [Burgers] ( 1969| l if we assume 
|Wi| Us and the absence of magnetic fields, the basic equa¬ 
tions of diffusion are 

m,rs - tUsTt 

Ms + m, 
(89) 


'"^Ps - Psg - Pes^ = ^ Kst(w, - W,) + ^ KstZst 


^ 2 5 _ 1 

= --Kssz'',rs - -^Ks,Zst- 
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ti^s 
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3ml + mjz'„ 
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+ -- 


- m, 

ntsfrit 


(iris + m,)^ 5 (Ms + Mt)^ 


msin, 


(ms + m,)2 


3 + Zst ^ I 


(90) 


where E is the quasi-static electric field and p^s is the average 
charge density of species s. These equations are still general, 
with the form of the resistance coefficients not yet fully speci¬ 
fied. The physics of the particular types of interactions within 
ideal gases is fully contained in the coefficients Kst, Zst, z'st, 
z's't, and z"/. 

For ionized gases, the resistance coefficients require eval¬ 
uation of collision integrals that diverge for a pure Coulomb 


potential. However, since the two-particle interaction poten¬ 
tial is only truly applicable on short length scales, an integra¬ 
tion cutoff or screened potential is commonly adopted. Burg¬ 
ers chooses to calculate resistance coefficients using a pure 
Coulomb potential truncated at the Debye radius 


Rd = 


ry2 


(91) 


which is assumed to be much larger than the inter-ion spac¬ 
ing. Indeed, for a plasma of one species, f^o/dion = (3r)“'^^. 
Applying this form of interaction to the collision integrals, 
the / - 1 integrals defined in Equation ( [84l l can be evalu¬ 
ated ( |Baalrud & Daligault|2014| l 


c(l) ^ . 

A2v4 
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(92) 


where A,, = yUj,a^j^/(ZjZ,e^). In order to perform the in¬ 
tegral in Equation ( |85l ). Burgers notes that the dependence of 
5® on V inside the logarithmic term is weak, so that we can 
replace there with its average value (v^) = hk^Tjpst- As¬ 
suming a very dilute plasma, so that » 1, Burg¬ 

ers then writes 


:(i) 




In 


hk^TRY) 


ZsZ,e^ 

and the final result for the resistance coefficients follows as 
\()^|nnsn,Z]Z}e^ 


(93) 
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(95) 


With these coefficients now fully specified. Burgers diffusion 
equations along with constraints such as charge neutrality and 
current neutrality form a closed set of equations, which can be 
solved for Wj, r^, E, and g from the input of a stellar profile. 

9.1.2. mesa’s Implementation of Thoul et al.’s Approach 

The diffusion routi ne originally imple mented in MESA was 
based on the work of Thoul et al. \\994\ . They start with the 
Burgers equations, writte n in a compact notation following 
Noerdlinger ( |1977 1978| l that is equivalent to Equatio ns ([89]l 
and (|90|l in one dimension. However, the approach of|Thoul| 
et al.|(|1994|l differs from Burgers’ original treatment in one 
important respect: the resistance coefficients are based on a 
modified result for the collision integrals. They follow Equa¬ 
tion ( |95] ) for the various Zst coefficients, which uses a pure 
Coulomb potential with a cutoff at the Debye length, but the 
Kst coefficients were derived from an alternative fitting of 
the Co ulomb logarithms introduced by |Iben & MacDonald 


(|1985|). Eor these coefficients, they define A - max(/?D, Tion), 
and use 


Ks, = 


16 nsn,Z]z}e^ 


Psta], 


1.6249 
X —-— In 


1-1-0.18769 


4kY2.TA 

ZsZ,e^ 


1.2 


(96) 


This expression is a fit to the numerical results of Eontaine & 
Michaud (19791, motivated by white dwarf conditions where 
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Burgers’ approximations for dealing with Equation ( |92] i are 
not valid (E > 1). Since this fit focuses on the strong coupling 
regime, and differs from Equation ( |94l i, these results can be 
incorrect in t he limit of a dilute plasma as we discuss later. 
Nevertheless, |Thoul et al.|(|1994|) elected to use Equation (|96ll 
under all conditions, since it provides an approximately cor¬ 
rect solution in a convenient closed form. 

Using Equations ( [89l l and ( |90l l along with the constraints of 
current neutralit y {Y^sPesWs - 0) an d local mass conservation 
(ZiPitVj = 0), Thoul et al. ( 1994[ ) express an entire closed 
system of equations in a dirnensionless matrix form suitable 
for numerical evaluation; 


p ( d\vip d\viT 

dr dr 


Z 

i=i 


d\nC, 


Jij- 


dr 


2S+2 


Z ^uWj, (97) 


J=i 


where S is the total number of species in the gas (includ¬ 
ing electrons) an d Cj = njjn^ is the concentration of the jth 
species. Consu lt j^oulet aTT (19941 for definitions of Kq, a,. 


Vi, jij, and Aij. The dehnition of W 


Wj^ 


Wj for j - 1.. .S, 

rj fovj-S + l. . .25, 

KQ^iiecE for j — 2S + 1, 
KQ^rigin^g for j - 2S + 2. 


(98) 


This is the vector containing the unknown quantities solved 
for after specifying Kq, a,, v,, y,^, and A,y. The routine pro¬ 
vided bylThoul et al.|(|1994|l inverts Equation (|97]l for one term 
in the left hand side at a time so as to find the “generalized 
diffusion coefficients,” which can be used to construct diffu¬ 
sion velocities or contributions from pressure, temperature, or 
concentrations individually. 

9.1.3. Modified Coefficients and Radiative Levitation as 
Implemented by Hu et al. 

|Hu et al.| ( |201 1[ ) extend the methods of |Thoul et al. ( |1994| l 
by introducing some key modifications. Eirst, they include 
an extra force term due to radiative levitation, so that Equa¬ 
tion ([89| becomes 


dps 

-T~ + - 5rad,i) - n.ZseE 

dr 

= Xi - tt's) + ^ KstZst 


m,rs - m,r, 


(99) 


t^S 


t^S 


OTo 


• m, 


where grad.s refers to the radiative acceleration on species s. 
Zs is the average charge of species s, allowing an account of 
partial ionization so that MjZje = p^s- They do not modify 
Equation (|90l)p^ 


In contrast to Thoul’s original ro utine, Hu et al. ( 2011 1 
use the resistance coefficients from [Paquette et al.| ( |1986| l, 
which were generated based on substantial improvements to 
Eontaine & Michaud ( 1979| l. In evaluating the collision inte- 
grals, [Paquette et al.|( ll)l56|l use a screened Coulomb potential 
of the form 

1/ r ^ - 7 7 2exp(-r/T) 

^st^P) — ZfZiG , ( 100 ) 


As written in equation (3) of[Hu et al.[(|2011^, their expression has two 
errors in the first term on the right hand side ot the hrst line: the sign is wrong, 
and it is missing resistance coefficients Kjj. Since neither of these en'ors 
propagates into later sections of the paper, it appears t hat both are simply 
typos, and otherwise their expression matches Equation j901 exactly. 


where, once again, A = max(/?D, dion). As we note below, this 
choice of A makes a substantial difference in strongly coupled 
plasmas, where the Debye radius no longer corresponds to a 
distance at which other nearby charged particles can signifi¬ 
cantly screen the Coulomb field. After se tting up the algebr a 
for a matrix solution very similar to that of|Thoul et al.[([1994[|, 
Hu et a l. (20TT) solve for the vector Wj (as dehned in Equa- 
tion[98]rappearing in the equation 

P I aiifiig^^ij d\np 

Ko\ k^T dr 


d\nT ^ d\nCj\ (101) 

;=i ' 7=1 


Many of the quan tities appear i ng in this equation are defi ned 
differently than in[Thoul et al.[([1994|l; see[Hu et al.[(|201 l|l for 
details. We can also solve this equation directly for the vector 


Wj to obtain 


^25 + 1 ^0 eE 


= -, ( 102 ) 


^ 25+2 K^^^tleMpg nipg’ 
the strength of the electric field relative to gravity. 

9.2. Analytic Expression for the Electric Field 

In some simple cases. Burgers equations can be solved to 
yield an analytic expression for the electric field, providing a 
usefu l test for MESA. Starting directly with his diffusion equa¬ 
tions, [Burgers] (|T969l) arrives at the following expressions for 
a pure plasma of electrons along with one species of ions 
(charge Ze): 
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(104) 


where w = Wi - We. Eor a plasma with only one ion species 
in diffusion equilibrium, the constraints of current neutrality 
and local mass conservation give w = 0. In the case of a pure 
hydrogen plasma, p - 2p^, and in hydrostatic equilibrium 
Vpe = V/7/2 = pg/2. Hence, we can solve the above set of 
equations to find 
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eE = --mpg - 


3 /2/:eeZ'; 


2'P*’ 2\5 Eie 


ksVT. 


(105) 


The coefficient for the temperature gradient term depends di¬ 
rectly on the nature of the resistance coefficients in the Burg¬ 
ers formalism, so different models of Coulomb collisions in 
ionized plasma will lead to different results for the electric 
field. 

As a slight generalization of Equation ( |105[ l in one dimen¬ 
sion, we write 


eE 1 ks dT 
Mpg 2 " Mpg dr 


(106) 


If we calculate the coefficient ag using the Burgers’ formalism 
with Equations (|95l) and (|9^, we find 


3/2 EeeZee , 

Ole^ -\- — - + Zi, 
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= 0.804 


(107) 
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A compa rable analytic expressio n for the electric field is pro¬ 
vided by |Roussel-Dupre| (|1981|, who applies a Boltzmann- 
Fokker-Planck approach to hnding diffusion coefficients for 
trace elements in hydrogen plasma. His treatment of diffu¬ 
sion is more precise than the Burgers’ formalism, but has the 
limitation of only being applicable in the case of nearly pure 
hydrogen with a diffusing trace element. His r esult for the 
electric field matches the form of Equation ( |106| l with the co¬ 
efficient Oe = 0.703. This provides another useful point of 
comparison in the specific case of nearly pure hydrogen plas¬ 
mas. Below we use this analytic expressio n as a test of th e 
updated resistance coefficients employed by Hu et al. (|2011[). 


9.3. Results and Comparisons 

We have constructed several simple MESA test cases in or¬ 
der to illustrate the effects of radiative levitation and differ¬ 
ent resistance coefficients. Where possible, we compare MESA 
output to corresponding analytic expressions. 


9.3.1. Electric Fields 

By default, MESA use s the resistance coefficients provided 
by Paquette et al. ( 1986 1 , but it can also use the resistance co- 
efhcients dehned by|lben & MacDonald|(|19851l, given here in 
Equation ( |96l l. In the case of a pure hydrogen star, the coeffi¬ 
cients given in Equation ( |96] ) lead directly to Equation ( |107| i, 
so these coefficients are especially useful in performing sim¬ 
ple comparisons of MESA output to a corresponding analytic 
expression. Due to the complicated numerical methods us ed 
to obtain the resistance coefficients ofjPaquette et al.|(| 198611, it 
is not possible to write down a directly corresponditig closed 
form analytic expression for the electric field, but results 
based on t hese more precise calc ulations compare favorably 
to those of |Rous sel-Dupre | (| 1981 |l in the case of a pure hydro¬ 
gen plasma. Starting with the MESA test suite, we constructed 
a solar mass pure hydrogen star, and we ran just long enough 
to turn on the diffusion routine and gather output for electric 
and gravitational fields. For such a star, we can compare MESA 
results for the elect ric fi eld directly to the analytic expression 
given in Equation ( |106[ l, with = 0.804 in the solutio n of 
Burgers (I1969|l and Og - 0.703 for Roussel-Dupre ( |1981| l. 

Figure p^l^ots the result of Equation ( |106| l fornoth values 
of a g, alo ng with the results from the diffusion routine (Equa¬ 
tion [TO^ for each type of resistance coefficients available in 
MESA. As expected, the curve calculated from the MES A dif¬ 
fusion routine output using the resistance coefficients of Iben 


& MacDonald 
with Ug - 0.8(lT as 
similar coefficients 


(1985 I closely matches the analytic expression 
IS calculated by |Burgers| (1969]) using hi 
. When using the more detailed numerics 


calculations fo r the resistance coefficients provided by Paque- 
|tte et al.|(1986| l, the diffusion routine output cl osely resembles 
the more precise analytic calculation given bylRoussel-Dupre 
( fT98T] ). 

The Sun provides another interesting test case for compar¬ 
ing the effects of using different resistance coefficients. An ex¬ 
ample solar model from the MESA test suite was run with dif¬ 
ferent choices of the resistance coefficients. Pigure|40] shows 
a slight difference bet ween the electric field s trengths relative 
to gravit y given by the Paquette et al. ( 1986[ l coefficients and 
those by Iben & MacDona 3 (|1985p. 


9.3.2. Gravitational Fields 

The MESA diffusion routine treats both the electric field and 
local gravitational acceleration as unknown quantities. MESA 



Figure 39. Compaiison of electric field strengths relative to gravity in a pure 
hydrogen star {M = 1.0M©, = 5.74 x 10^ K, L = 0.576Lo) with nu¬ 

clear burning artificially suppressed in the MESA routine to avoid any helium 
cont amina tion. Solid lines represent the analytic expression given by Equa¬ 
tion for two different values of the coefficient ag. Dashed lines repre¬ 
sent output from the MESA diffusion routine as described in Equation jl02) , 
with the only difference being the resistance coefficients used to soWe the 
Burgers equations. 
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Figure 40. Comparison of electric field strengths relative to gravity using 
different resistance coefficients in a solar model. 


records the quantity W 2 S +2 (Equation [98]l, used to calculate 
the gravitational acceleration from the dltfusion routine; 


.gdiff 


KqW2S+2 

n^m^ 


(108) 


This expression for gditt is independent of the simpler expres¬ 
sion for local gravitational acceleration ggauss = Gmfr^. Eig- 
ure 41 compares ggauss and gdift for a typical profile found 


using the example solar model from the MESA test suite. In 
Eigure 41a profile from a star of larger mass (M =1.5 Mq) 
shows disagreement between the gravity outputs in the con¬ 
vective core. The Burgers formalism assumes heat transfer 
that is correlated with temperature gradients through the resid¬ 
ual heat flow vectors defined in Equation ( |88] l. This assump¬ 
tion breaks down when most of energy is transported by con¬ 
vection; however, the effects of diffusion in this region are 
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completely overwhelmed by convective mixing and are there¬ 
fore inconsequential. 



Figure 41. Comparison of gravitational fields obtained from gag and 
ggauss in two MESA test suite cases. The two lines representing the Sun 
(age = 4.57 Gyr) show good agreement, while the two lines representing a 
1.5 Mq star disagree in regions with large convective flux where diffusion is 
inconsequential. 


9.3.3. Radiative Levitation 
mesa’s implementation of radiative levitation is based on 


Hu et al.| (|2011||. Figure |42| shows an abundance profile of 
a subdwarr' B star model produced by MESA, where radiative 
levitation is responsible for the presence of ^®Fe, ^^Ni, and 


other metals near the surface (as also seen in figure 3 of Hu 
let al.|2011| l. 


9.3.4. White Dwarf Sedimentation 

In a cooling WD, diffusion governs sedimentation over long 
timescales. The assumptions behind the formalism of the 
Burgers equations do not hold under white dwarf conditions: 

• The Burgers equations assume all particle species sat¬ 
isfy an ideal gas equation of state. In the context of a 
degenerate WD both electrons and ions violate this as¬ 
sumption. 

• The very dense, strongly coupled (F > 1) conditions of 
a WD call into question the validity of the two-particle 
scattering picture used to calculate the ion resistance 
coefficients. 


Nevertheless, for lack of a better option, previous studi es have 
relied on the Burgers equations with the coefficients of Paque- 


|tte et al.| (|1986|). For example, see Corsico et al. (2002 1 . 

Figure |4^¥ hows an abundance prohle produced by MESA 
for a CO WD after 4 Gyr of evolution, where diffusion gov¬ 
erns sedimentation in the outer layers. The vertical lines in 
Figure [4^mark the outer boundaries of regions where the two 
concernslisted above become significant. Nearly all of the 
WD resides inside at least one of these regimes, and much 



Figure 42. Abundance profile of a subdwarf B star model (M = 0.462 Mq, 
Teff = 2.67 X 10"^ K, L = 1.12Lq, age = 5 Myr) showing the etfects of radia¬ 
tive levitation with a layer of ^®Fe/^®Ni at the surface. 


of the interesting diffusion sedimentation occurs inside re¬ 
gions that are both significantly coupled and highly degen¬ 
erate. Thus, improvements to the treatment of diffusion are 
clearly necessary before we are able to describe diffusion in 
WDs adequately. This MESA run turns off diffusion for F > 50, 
where we expect strong coupling to substantially modify the 
underlying equations. 




Figure 43. Abundance profile of a CO WD {M = 0.611 M©, Feff = 5.16 x 
10'^ K, L = 9.29 X 10“^ Lq) after 4 Gyr of WD evolution. The region left of 
the blue, dashed line is the interior of the WD, where T > 1. Left of the red, 
dashed line T > 50, and diffusion has been turned off for this region. The 
electrons are an ideal gas to the right of the black dot-dash line. 


9.4. Expanding the Domain of Validity and Next Steps 
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The validity of the Boltzmann approach be comes question 
able as r > 1 and the ions become a liquid. |Bildsten & Hall 
(20011 estimated the diffusion coefficient in this liquid regime 
by using the Stokes-Einstein relation. However, for a broad- 
based code such as MESA, we need to implement diffusion into 
the r > 1 regime in a manner that allows for a smooth transi- 
ti on between cou pling regimes. 

Paquette et al. (|1986|l successfully described diffusion in a 


regime of intermediate coupling through the use of screened 
potentials, which are a way to account for the collective nature 
of interactions in a dense plasma. Though there is no rigorous 
reason to expect that a formalism based on the two-particle 
scattering picture should work well as T —» 1, their compari¬ 
son to simulations verified that this description of diffusion is 
very accurate for T < 1. 

Can these approximatio ns be extrapolated to the stro ngly 
coupled regime of T > 1 ? |Baalrud & Daligault| (|2013|l pro¬ 
vide a method for numerically calculating resistance coeffi¬ 
cients using a hypernetted chain (HNC) approximation from 
effective potentials. Figure l44| compares their HNC results 
(diamonds) to their MoleculaFDynamics (MD) simulations of 
a one-component plasma (OCP, circles) for the self-diffusion 
coefficient D*, defined by 


D* = 


D 


atr, 
ion P 


(109) 


where (Op is the plasma frequency and D - 20^^ (the factor 
of 2 in this definition ensures that if we redefine species s 
in terms of two subspecies si and S 2 , then D = The 

general expression for the interdiffusion coefficient is 




tlstlt 


k^T 


iis + n, Ksiil - A)’ 


( 110 ) 


where the 1 - A term in the denominator accounts for a second 
order correction that can be defined using 


A = 






55r 


( 11 ) 


• 2or 


,( 12 ) 




( 111 ) 


For reference, we also include a direct fit of Daligault & 
Murillo|(|2005|) to the MD data of|Ranganathan k al. (|i003|l, 
given by 
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D* = 0.0028 H- 0.005251^-1 


1.154 


( 112 ) 


The agreement between the HNC and MD simulations shows 
that the HNC does a better job of accounting for correlation 
physics in strongly coupled plasmas than a simple screened 
Coulomb potential and allows for a surprising (and still phys¬ 
ically unexplained) extension of the Burgers formalism into 
the strongly coupled regime. This recent work allows us to 
go into the large F limit with the Burgers formalism, but the 
question remains as to how we obtain diffusion coefficients in 
a reliable manner. 

The self-diffusion coefficients from the two options in MESA 
are shown in Figure 44 and correlate with the MD data better 
than expected for thenigh F regime. In particular, the agree- 
ment is much better than that shown in figure 2 ofjBaalrud & 
Daligault] (|2013|l for either “cutoff” or “screened” Coulomb 
methods. The reason for this agreement is that both MESA 
implementations use the inter-ion spacing rather than the De¬ 
bye length once F > 1/3, which yields favorable scalings in 



Figure 44. Compilation of the self-diffusion coefficients obtaine d from dif- 
ferent methods. “MD Data” and “HNC” points are taken from|Baalrud &| 
IDaligaultlpOlSj. The s olid black line is the result of the MESA calculation 
using the coefficients of|Paquette et al. 1(1986). The dashed green line is the 
result of the calculation us ing the resistance coefficients from the original rou- 
tine ofIT houl et alT|l l 994| based on the fit to the Coulomb logarithm found in 
jlben & Macnonffl|(iy!15i, given here in Equation {96\. The dashed purple 
line represents the ht to MD data given here in Equation (112). 


the high F limit. Iben & MacDonald (1985 i constructed their 


fitting formula based on a few numerical results for F > 1. IPa- 
quette et al. ( 1986]) also showed that their formalism can be 


extended to F > 1 as long as the inter-ion spacing is used 
rather than the Debye radius for the screening length. 

Though MESA does not yet provide the capability of imple¬ 
menting resistance coefficients based on the HNC method, we 
hope to accomplish this in the near future by means of a table 
similar to that provided for the coefficients ofjPaquette et al.j 
(|1986j). For a more thorough discussion of these methods and 
the likely path of a pplication to mixtures, consult |Beznogov| 
|& Yakovlev| ( |2014) . We will also need to correctly account for 
the electron degeneracy and the non-ideal equation of state for 
the ions, both of which modify the electrostatic field needed 
to correctly determine the forces that drive diffusion. 


10. SOFTWARE INFRASTRUCTURE 

Here we describe a number of changes to MESA that have 
occurred since Paper II and are of potential interest to users of 
MESA or developers of similar software. 

MESA can be compiled with either the GNU or Intel For¬ 
tran compilers, runs on multiple operating systems (Windows, 
OS X, and Finux), and can use different numbers of OpenMP 
threads. It is necessary to regularly test that the code is per¬ 
forming correctly across the different combinations of com¬ 
piler, OS, thread count. To this end, developers and engaged 
users run the MESA test suite on a wide range of systems before 
each release. 

Previously, test cases in the MESA test suite accepted dif¬ 
ferent results so long as they were within a certain tolerance, 
an appropriate choice for testing scientific results where the 
physical uncertainties are much greater than the numerical 
ones. However, we found that this made detecting and track¬ 
ing bugs across platforms difficult. For the purposes of code 
testing, it is much better to insist that any inconsistency is a 
problem, no matter how small. 

Motivated by this challenge, MESA now provides bit-for-bit 
consistency for all results across all the supported platforms. 
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It is essential to emphasize that the goal of this achievement 
is to enable better testing. It allows users to exactly repro¬ 
duce the results of others, independent of platform differ¬ 
ences, which is especially useful to developers attempting to 
reproduce bugs. The achievement of bit-for-bit consistency is 
not a claim that the results of MESA calculations are physically 
accurate or numerically converged to any specific degree. 

This bit-for-bit consistency was achieved via the following 
choices: 

• Using parallel algorithms that give identical results in¬ 
dependent of number of threads or order of thread 
execution. MESA’s linear algebra solver is based on 
BCYCLIC ( jHirshman et al.||2010| l. It sub-divides the 
work between threads based on the the size of the ma¬ 
trix rather than on the number of threads available. It 
is also necessary to avoid OpenMP reduction clauses, 
which provide no guarantees on ordering of operations. 

• Specifying compiler flags that forbid the compiler from 
making any optimization that can affect floating point 
precision (e.g., forbid re-association and fast math op¬ 
erations). Most optimizations are still allowed. 

• Using an I/O library that does precise conversion from 
binary to ASCII for double precision numbers. 

• Using a math library that gives consistent results for 
operations such as log, exp, sin, cos, pow. MESA uses 
CRLIB1VC3 in round towards zero mode. The choice 
to use a math library that gives exact results is not be¬ 
cause 16 digit accuracy from the math routines in nec¬ 
essary. Rather, we want consistent results across sup¬ 
ported platforms and this is the best way to achieve this 
consistency. 

• Replacing integer power expressions (i.e., x**3) by re¬ 
peated multiplications (i.e., x*x*x). Different compil¬ 
ers implement integer powers differently, giving differ¬ 
ent results. 

Having achieved bit-for-bit identical results, we can test 
files for exact equality. This applies both to the module-by- 
module tests that run at installation time and the case by case 
tests in the star and binary test suites. These test cases com¬ 
pare the final model from the test run to a saved result from a 
previous MESA version. If they are not exactly the same, the 
test fails. The test is also restarted from an intermediate state 
to confirm that runs which are stopped and restarted yield ex¬ 
actly the same results as those that are not. 

While MESAstar is parallelized via OpenMP, the install 
process has historically been serial. MESA contains approxi¬ 
mately 1000 Fortran files and so the ability to compile more 
than one file simultaneously has the potential to provide sig¬ 
nificant reductions in the time needed to install MESA. Re¬ 
cently the compilation step has been parallelized, enabled by 
the automated dependency generation tool makedep£90{^al- 
lowing multiple instances of the Fortran compiler to be in¬ 
voked simultaneously. This is of particular utility for devel¬ 
opers who may recompile MESA frequently. 

Since Paper II the main MESA websitq^has undergone sig¬ 
nificant revision, making it easier for new users to get started 

http://lipforge.ens-lyon.fr/www/crlibm/index.html 

http://personal.inet.fi/private/erikedelmann/makedepf90/ 

** http://mesa.sourceforge.net 


with MESA. This restructuring has also made it easier for the 
developers to keep material up-to-date as MESA evolves. One 
of the most important improvements is that the files that doc¬ 
ument the default value of each MESA option use the Mark- 
dowrp^ markup language. This allows documentation web 
pages to be generated automatically for each MESA release. 

Improvements have also been made to the distribution of 
MESA. Previously, MESA was available only by checking out 
the source code using the Subversiorj^ version control sys¬ 
tem. Now, every release version of MESA (including past re¬ 
leases) is available for download as a ZIP archive. This is 
simpler and saves bandwidth and disk space. It has quickly 
become the preferred way to install MESA with the ZIP file of 
the current release being downloaded tens of times per week. 


11. SUMMARY AND CONCLUSIONS 

We have explained and, where possible, verified or vali¬ 
dated, major new capabilities and improvements implemented 
in MESA since the publication of Paper I and Paper II. These 
advancements include interacting binary systems (Section]^, 
implicit hydrodynamics and shocks (Section ffl, in-situ usage 
of large reaction networks, especially for X^ray bursts and 
core-collapse supernova progenitors (Section [^, and the ex¬ 
plosion of massive stars (Section]^. These new capabilities 
will allow for extended exploration of core collapse progen¬ 
itors and the sensitivity of shock nucleosynthesis to their ex¬ 
plosion mechanism. The full coupling of MESA to the GYRE 
non-adiabatic pulsation instrument (Section!^ has already re¬ 
vealed the richness of the instability strips ror massive stars 
and enables the continued growth of astero-seismology across 
the HR diagram. Progress in the treatment of mass accre¬ 
tion (Section]^ and weak reaction rates (Sectionwill im¬ 
prove studies of their impact on stellar evolution. We also dis¬ 
cuss the domain of validity for particle diffusion within MESA 
and describe a path forward for extending diffusion into the 
regime relevant to WD sedimentation (Section]^. We also de¬ 
scribe significant improvements to the infrastructure of MESA 
(Section [T0|). MESAstar input files and related materials for 
all the figures are available at http://mesastar.org 

These hitherto unpublished advancements have already en- 


abled a number of studies in interacting binary systems ( 

Wolf 

|et al.|2013[|Pavlovskii & Ivanova|2015[|Vos et al.|2015 

1 and 

stellar pulsations (jPapics et al.|2014 

jStello et al.|2014 Quinn 

et al. 

2015[|Cunha et al. 2015 1 , an 

d led to the discovery of 


ONeMg cores towards AIC ( Schwab et al.|20f5] l. It also en¬ 
abled the first three dimensional simulations of the final min¬ 
utes of iron core growth in a massive star up to and including 


the point of core gravitational instability and collapse ( Couch] 
et al.|201? I. In addition, these enhanced capabilities have al¬ 


lowed for applications of MESAstar that were not initially 
envisioned, such as the treatment of M agneto-Rotational In¬ 
stability in stars ( Wheeler et al.||20T5]l, effects of axions on 
nucleosynthesis in massive stars ( Aoyama & SuzuE||2015|), 
and particle physics beyond the Standard Model (|Curtin & 
|Tsai|2014| ). 

As a community software instrument for stellar astro¬ 
physics new directions for MESA will be driven by: features 
useful to the MESA user community, advances in the physics 
modules, algorithmic developments, and architectural evolu¬ 
tion. Potential examples include a treatment of ionization in 


http://daringfireball.net/projects/markdown/ 

https://subversion.apache.org/ 





































42 


the equation of state for an arbitrary composition across an ex¬ 
panded region in the p-T plane, non-linear pulsations, Monte 
Carlo based thermonuclear reaction rates, modules for sub¬ 
sonic flame propagation, ports to additional architectures, and 
a web-interface to MESA for education. 
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